/*
    Uncmin_f77.java copyright claim:

    This software is based on the public domain UNCMIN routines.
    It was translated from FORTRAN to Java by a US government employee 
    on official time.  Thus this software is also in the public domain.


    The translator's mail address is:

    Steve Verrill 
    USDA Forest Products Laboratory
    1 Gifford Pinchot Drive
    Madison, Wisconsin
    53705


    The translator's e-mail address is:

    steve@ws13.fpl.fs.fed.us


 ***********************************************************************

DISCLAIMER OF WARRANTIES:

THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND. 
THE TRANSLATOR DOES NOT WARRANT, GUARANTEE OR MAKE ANY REPRESENTATIONS 
REGARDING THE SOFTWARE OR DOCUMENTATION IN TERMS OF THEIR CORRECTNESS, 
RELIABILITY, CURRENCY, OR OTHERWISE. THE ENTIRE RISK AS TO 
THE RESULTS AND PERFORMANCE OF THE SOFTWARE IS ASSUMED BY YOU. 
IN NO CASE WILL ANY PARTY INVOLVED WITH THE CREATION OR DISTRIBUTION 
OF THE SOFTWARE BE LIABLE FOR ANY DAMAGE THAT MAY RESULT FROM THE USE 
OF THIS SOFTWARE.

Sorry about that.

 ***********************************************************************


History:

Date        Translator        Changes

4/14/98     Steve Verrill     Translated

 */

package doglegMethod;

/**
 * 
 *<p>
 * This class contains Java translations of the UNCMIN unconstrained
 * optimization routines. See R.B. Schnabel, J.E. Koontz, and B.E. Weiss, <i>A
 * Modular System of Algorithms for Unconstrained Minimization</i>, Report
 * CU-CS-240-82, Comp. Sci. Dept., University of Colorado at Boulder, 1982.
 * 
 *<p>
 * <b>IMPORTANT:</b> The "_f77" suffixes indicate that these routines use
 * FORTRAN style indexing. For example, you will see
 * 
 * <pre>
 *   for (i = 1; i <= n; i++)
 *</pre>
 * 
 * rather than
 * 
 * <pre>
 *   for (i = 0; i < n; i++)
 *</pre>
 * 
 * To use the "_f77" routines you will have to declare your vectors and matrices
 * to be one element larger (e.g., v[101] rather than v[100], and a[101][101]
 * rather than a[100][100]), and you will have to fill elements 1 through n
 * rather than elements 0 through n - 1. Versions of these programs that use
 * C/Java style indexing will eventually be available. They will end with the
 * suffix "_j".
 * 
 *<p>
 * This class was translated by a statistician from a FORTRAN version of UNCMIN.
 * It is NOT an official translation. It wastes memory by failing to use the
 * first elements of vectors. When public domain Java optimization routines
 * become available from the people who produced UNCMIN, then <b>THE CODE
 * PRODUCED BY THE NUMERICAL ANALYSTS SHOULD BE USED</b>.
 * 
 *<p>
 * Meanwhile, if you have suggestions for improving this code, please contact
 * Steve Verrill at steve@ws13.fpl.fs.fed.us.
 * 
 *@author (translator)Steve Verrill
 *@version .5 --- April 14, 1998
 * 
 */

@SuppressWarnings("all")
public class Uncmin_f77 extends Object {

	/**
	 * 
	 *<p>
	 * The optif0_f77 method minimizes a smooth nonlinear function of n
	 * variables. A method that computes the function value at any point must be
	 * supplied. (See Uncmin_methods.java and UncminTest.java.) Derivative
	 * values are not required. The optif0_f77 method provides the simplest user
	 * access to the UNCMIN minimization routines. Without a recompile, the user
	 * has no control over options. For details, see the Schnabel et al
	 * reference and the comments in the code.
	 * 
	 *Translated by Steve Verrill, August 4, 1998.
	 * 
	 *@param n
	 *            The number of arguments of the function to minimize
	 *@param x
	 *            The initial estimate of the minimum point
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param xpls
	 *            The final estimate of the minimum point
	 *@param fpls
	 *            The value of f_to_minimize at xpls
	 *@param gpls
	 *            The gradient at the local minimum xpls
	 *@param itrmcd
	 *            Termination code ITRMCD = 0: Optimal solution found ITRMCD =
	 *            1: Terminated with gradient small, xpls is probably optimal
	 *            ITRMCD = 2: Terminated with stepsize small, xpls is probably
	 *            optimal ITRMCD = 3: Lower point cannot be found, xpls is
	 *            probably optimal ITRMCD = 4: Iteration limit (150) exceeded
	 *            ITRMCD = 5: Too many large steps, function may be unbounded
	 *@param a
	 *            Workspace for the Hessian (or its estimate) and its Cholesky
	 *            decomposition
	 *@param udiag
	 *            Workspace for the diagonal of the Hessian
	 * 
	 * 
	 */

	public static void optif0_f77(int n, double x[], Uncmin_methods minclass,
			double xpls[], double fpls[], double gpls[], int itrmcd[],
			double a[][], double udiag[]) {

		/*
		 * 
		 * Here is a copy of the optif0 FORTRAN documentation:
		 * 
		 * SUBROUTINE OPTIF0(NR,N,X,FCN,XPLS,FPLS,GPLS,ITRMCD,A,WRK) c implicit
		 * double precision (a-h,o-z) c C C PURPOSE C ------- C PROVIDE SIMPLEST
		 * INTERFACE TO MINIMIZATION PACKAGE. C USER HAS NO CONTROL OVER
		 * OPTIONS.
		 * 
		 * 
		 * C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C X(N) --> INITIAL ESTIMATE OF MINIMUM C FCN -->
		 * NAME OF ROUTINE TO EVALUATE MINIMIZATION FUNCTION. C MUST BE DECLARED
		 * EXTERNAL IN CALLING ROUTINE. C XPLS(N) <-- LOCAL MINIMUM C FPLS <--
		 * FUNCTION VALUE AT LOCAL MINIMUM XPLS C GPLS(N) <-- GRADIENT AT LOCAL
		 * MINIMUM XPLS C ITRMCD <-- TERMINATION CODE C A(N,N) --> WORKSPACE C
		 * WRK(N,9) --> WORKSPACE C
		 */

		int msg[] = new int[2];

		double typsiz[] = new double[n + 1];
		double g[] = new double[n + 1];
		double p[] = new double[n + 1];
		double sx[] = new double[n + 1];
		double wrk0[] = new double[n + 1];
		double wrk1[] = new double[n + 1];
		double wrk2[] = new double[n + 1];
		double wrk3[] = new double[n + 1];

		double dlt[] = new double[2];
		double fscale[] = new double[2];
		double stepmx[] = new double[2];

		int ndigit[] = new int[2];

		int i, ig, it;
		int lt;

		int method[] = new int[2];
		int iexp[] = new int[2];
		int itnlim[] = new int[2];
		int iagflg[] = new int[2];
		int iahflg[] = new int[2];

		double gradtl[] = new double[2];
		double steptl[] = new double[2];

		Uncmin_f77.dfault_f77(n, x, typsiz, fscale, method, iexp, msg, ndigit,
				itnlim, iagflg, iahflg, dlt, gradtl, stepmx, steptl);

		Uncmin_f77.optdrv_f77(n, x, minclass, typsiz, fscale, method, iexp,
				msg, ndigit, itnlim, iagflg, iahflg, dlt, gradtl, stepmx,
				steptl, xpls, fpls, gpls, itrmcd, a, udiag, g, p, sx, wrk0,
				wrk1, wrk2, wrk3);

		if (itrmcd[1] == 1) {

			System.out.print("\nUncmin WARNING --- itrmcd = 1, "
					+ "probably converged, gradient small\n\n");

		} else if (itrmcd[1] == 2) {

			System.out.print("\nUncmin WARNING --- itrmcd = 2, "
					+ "probably converged, stepsize small\n\n");

		} else if (itrmcd[1] == 3) {

			System.out.print("\nUncmin WARNING --- itrmcd = 3, "
					+ "cannot find lower point\n\n");

		} else if (itrmcd[1] == 4) {

			System.out.print("\nUncmin WARNING --- itrmcd = 4, "
					+ "too many iterations\n\n");

		} else if (itrmcd[1] == 5) {

			System.out.print("\nUncmin WARNING --- itrmcd = 5, "
					+ "too many large steps, " + "possibly unbounded\n\n");

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The optif9_f77 method minimizes a smooth nonlinear function of n
	 * variables. A method that computes the function value at any point must be
	 * supplied. (See Uncmin_methods.java and UncminTest.java.) Derivative
	 * values are not required. The optif9 method provides complete user access
	 * to the UNCMIN minimization routines. The user has full control over
	 * options. For details, see the Schnabel et al reference and the comments
	 * in the code.
	 * 
	 *Translated by Steve Verrill, August 4, 1998.
	 * 
	 *@param n
	 *            The number of arguments of the function to minimize
	 *@param x
	 *            The initial estimate of the minimum point
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param typsiz
	 *            Typical size for each component of x
	 *@param fscale
	 *            Estimate of the scale of the objective function
	 *@param method
	 *            Algorithm to use to solve the minimization problem = 1 line
	 *            search = 2 double dogleg = 3 More-Hebdon
	 *@param iexp
	 *            = 1 if the optimization function f_to_minimize is expensive to
	 *            evaluate, = 0 otherwise. If iexp = 1, then the Hessian will be
	 *            evaluated by secant update rather than analytically or by
	 *            finite differences.
	 *@param msg
	 *            Message to inhibit certain automatic checks and output
	 *@param ndigit
	 *            Number of good digits in the minimization function
	 *@param itnlim
	 *            Maximum number of allowable iterations
	 *@param iagflg
	 *            = 0 if an analytic gradient is not supplied
	 *@param iahflg
	 *            = 0 if an analytic Hessian is not supplied
	 *@param dlt
	 *            Trust region radius
	 *@param gradtl
	 *            Tolerance at which the gradient is considered close enough to
	 *            zero to terminate the algorithm
	 *@param stepmx
	 *            Maximum allowable step size
	 *@param steptl
	 *            Relative step size at which successive iterates are considered
	 *            close enough to terminate the algorithm
	 *@param xpls
	 *            The final estimate of the minimum point
	 *@param fpls
	 *            The value of f_to_minimize at xpls
	 *@param gpls
	 *            The gradient at the local minimum xpls
	 *@param itrmcd
	 *            Termination code ITRMCD = 0: Optimal solution found ITRMCD =
	 *            1: Terminated with gradient small, X is probably optimal
	 *            ITRMCD = 2: Terminated with stepsize small, X is probably
	 *            optimal ITRMCD = 3: Lower point cannot be found, X is probably
	 *            optimal ITRMCD = 4: Iteration limit (150) exceeded ITRMCD = 5:
	 *            Too many large steps, function may be unbounded
	 *@param a
	 *            Workspace for the Hessian (or its estimate) and its Cholesky
	 *            decomposition
	 *@param udiag
	 *            Workspace for the diagonal of the Hessian
	 * 
	 */

	public static void optif9_f77(int n, double x[], Uncmin_methods minclass,
			double typsiz[], double fscale[], int method[], int iexp[],
			int msg[], int ndigit[], int itnlim[], int iagflg[], int iahflg[],
			double dlt[], double gradtl[], double stepmx[], double steptl[],
			double xpls[], double fpls[], double gpls[], int itrmcd[],
			double a[][], double udiag[]) {

		/*
		 * 
		 * Here is a copy of the optif9 FORTRAN documentation:
		 * 
		 * SUBROUTINE OPTIF9(NR,N,X,FCN,D1FCN,D2FCN,TYPSIZ,FSCALE, +
		 * METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR, +
		 * DLT,GRADTL,STEPMX,STEPTL, + XPLS,FPLS,GPLS,ITRMCD,A,WRK) c implicit
		 * double precision (a-h,o-z) c C C PURPOSE C ------- C PROVIDE COMPLETE
		 * INTERFACE TO MINIMIZATION PACKAGE. C USER HAS FULL CONTROL OVER
		 * OPTIONS. C C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX
		 * C N --> DIMENSION OF PROBLEM C X(N) --> ON ENTRY: ESTIMATE TO A ROOT
		 * OF FCN C FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION
		 * C MUST BE DECLARED EXTERNAL IN CALLING ROUTINE C FCN: R(N) --> R(1) C
		 * D1FCN --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE GRADIENT C OF
		 * FCN. MUST BE DECLARED EXTERNAL IN CALLING ROUTINE C D2FCN -->
		 * (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE HESSIAN OF C OF FCN. MUST
		 * BE DECLARED EXTERNAL IN CALLING ROUTINE C TYPSIZ(N) --> TYPICAL SIZE
		 * FOR EACH COMPONENT OF X C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE
		 * FUNCTION C METHOD --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
		 * C =1 LINE SEARCH C =2 DOUBLE DOGLEG C =3 MORE-HEBDON C IEXP --> =1 IF
		 * OPTIMIZATION FUNCTION FCN IS EXPENSIVE TO C EVALUATE, =0 OTHERWISE.
		 * IF SET THEN HESSIAN WILL C BE EVALUATED BY SECANT UPDATE INSTEAD OF C
		 * ANALYTICALLY OR BY FINITE DIFFERENCES C MSG <--> ON INPUT: (.GT.0)
		 * MESSAGE TO INHIBIT CERTAIN C AUTOMATIC CHECKS C ON OUTPUT: (.LT.0)
		 * ERROR CODE; =0 NO ERROR C NDIGIT --> NUMBER OF GOOD DIGITS IN
		 * OPTIMIZATION FUNCTION FCN C ITNLIM --> MAXIMUM NUMBER OF ALLOWABLE
		 * ITERATIONS C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED C IAHFLG -->
		 * =1 IF ANALYTIC HESSIAN SUPPLIED C IPR --> DEVICE TO WHICH TO SEND
		 * OUTPUT C DLT --> TRUST REGION RADIUS C GRADTL --> TOLERANCE AT WHICH
		 * GRADIENT CONSIDERED CLOSE C ENOUGH TO ZERO TO TERMINATE ALGORITHM C
		 * STEPMX --> MAXIMUM ALLOWABLE STEP SIZE C STEPTL --> RELATIVE STEP
		 * SIZE AT WHICH SUCCESSIVE ITERATES C CONSIDERED CLOSE ENOUGH TO
		 * TERMINATE ALGORITHM C XPLS(N) <--> ON EXIT: XPLS IS LOCAL MINIMUM C
		 * FPLS <--> ON EXIT: FUNCTION VALUE AT SOLUTION, XPLS C GPLS(N) <--> ON
		 * EXIT: GRADIENT AT SOLUTION XPLS C ITRMCD <-- TERMINATION CODE C
		 * A(N,N) --> WORKSPACE FOR HESSIAN (OR ESTIMATE) C AND ITS CHOLESKY
		 * DECOMPOSITION C WRK(N,8) --> WORKSPACE C
		 */

		double g[] = new double[n + 1];
		double p[] = new double[n + 1];
		double sx[] = new double[n + 1];
		double wrk0[] = new double[n + 1];
		double wrk1[] = new double[n + 1];
		double wrk2[] = new double[n + 1];
		double wrk3[] = new double[n + 1];

		// MINIMIZE FUNCTION

		Uncmin_f77.optdrv_f77(n, x, minclass, typsiz, fscale, method, iexp,
				msg, ndigit, itnlim, iagflg, iahflg, dlt, gradtl, stepmx,
				steptl, xpls, fpls, gpls, itrmcd, a, udiag, g, p, sx, wrk0,
				wrk1, wrk2, wrk3);

		if (itrmcd[1] == 1) {

			System.out.print("\nUncmin WARNING --- itrmcd = 1, "
					+ "probably converged, gradient small\n\n");

		} else if (itrmcd[1] == 2) {

			System.out.print("\nUncmin WARNING --- itrmcd = 2, "
					+ "probably converged, stepsize small\n\n");

		} else if (itrmcd[1] == 3) {

			System.out.print("\nUncmin WARNING --- itrmcd = 3, "
					+ "cannot find lower point\n\n");

		} else if (itrmcd[1] == 4) {

			System.out.print("\nUncmin WARNING --- itrmcd = 4, "
					+ "too many iterations\n\n");

		} else if (itrmcd[1] == 5) {

			System.out.print("\nUncmin WARNING --- itrmcd = 5, "
					+ "too many large steps, " + "possibly unbounded\n\n");

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The bakslv_f77 method solves Ax = b where A is an upper triangular
	 * matrix. Note that A is input as a lower triangular matrix and this method
	 * takes its transpose implicitly.
	 * 
	 *Translated by Steve Verrill, April 14, 1998.
	 * 
	 *@param n
	 *            Dimension of the problem
	 *@param a
	 *            n by n lower triangular matrix (preserved)
	 *@param x
	 *            The solution vector
	 *@param b
	 *            The right-hand side vector
	 * 
	 */

	public static void bakslv_f77(int n, double a[][], double x[], double b[]) {

		/*
		 * 
		 * Here is a copy of the bakslv FORTRAN documentation:
		 * 
		 * SUBROUTINE BAKSLV(NR,N,A,X,B)
		 * 
		 * C C PURPOSE C ------- C SOLVE AX=B WHERE A IS UPPER TRIANGULAR
		 * MATRIX. C NOTE THAT A IS INPUT AS A LOWER TRIANGULAR MATRIX AND C
		 * THAT THIS ROUTINE TAKES ITS TRANSPOSE IMPLICITLY. C C PARAMETERS C
		 * ---------- C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF
		 * PROBLEM C A(N,N) --> LOWER TRIANGULAR MATRIX (PRESERVED) C X(N) <--
		 * SOLUTION VECTOR C B(N) --> RIGHT-HAND SIDE VECTOR C C NOTE C ---- C
		 * IF B IS NO LONGER REQUIRED BY CALLING ROUTINE, C THEN VECTORS B AND X
		 * MAY SHARE THE SAME STORAGE. C
		 */

		int i, ip1, j;
		double sum;

		// SOLVE (L-TRANSPOSE)X=B. (BACK SOLVE)

		i = n;
		x[i] = b[i] / a[i][i];

		while (i > 1) {

			ip1 = i;
			i--;

			sum = 0.0;

			for (j = ip1; j <= n; j++) {

				sum += a[j][i] * x[j];

			}

			x[i] = (b[i] - sum) / a[i][i];

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The chlhsn_f77 method finds "THE L(L-TRANSPOSE) [WRITTEN LL+]
	 * DECOMPOSITION OF THE PERTURBED MODEL HESSIAN MATRIX A+MU*I(WHERE MU\0 AND
	 * I IS THE IDENTITY MATRIX) WHICH IS SAFELY POSITIVE DEFINITE. IF A IS
	 * SAFELY POSITIVE DEFINITE UPON ENTRY, THEN MU=0."
	 * 
	 *Translated by Steve Verrill, April 14, 1998.
	 * 
	 *@param n
	 *            Dimension of the problem
	 *@param a
	 *            On entry: A is the model Hessian (only the lower triangle and
	 *            diagonal stored) On exit: A contains L of the LL+
	 *            decomposition of the perturbed model Hessian in the lower
	 *            triangle and diagonal, and contains the Hessian in the upper
	 *            triangle and udiag
	 *@param epsm
	 *            Machine epsilon
	 *@param sx
	 *            Scaling vector for x
	 *@param udiag
	 *            On exit: Contains the diagonal of the Hessian
	 * 
	 */

	public static void chlhsn_f77(int n, double a[][], double epsm,
			double sx[], double udiag[]) {

		/*
		 * 
		 * Here is a copy of the chlhsn FORTRAN documentation:
		 * 
		 * SUBROUTINE CHLHSN(NR,N,A,EPSM,SX,UDIAG)
		 * 
		 * C C PURPOSE C ------- C FIND THE L(L-TRANSPOSE) [WRITTEN LL+]
		 * DECOMPOSITION OF THE PERTURBED C MODEL HESSIAN MATRIX A+MU*I(WHERE
		 * MU\0 AND I IS THE IDENTITY MATRIX) C WHICH IS SAFELY POSITIVE
		 * DEFINITE. IF A IS SAFELY POSITIVE DEFINITE C UPON ENTRY, THEN MU=0. C
		 * C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C A(N,N) <--> ON ENTRY; "A" IS MODEL HESSIAN
		 * (ONLY LOWER C TRIANGULAR PART AND DIAGONAL STORED) C ON EXIT: A
		 * CONTAINS L OF LL+ DECOMPOSITION OF C PERTURBED MODEL HESSIAN IN LOWER
		 * TRIANGULAR C PART AND DIAGONAL AND CONTAINS HESSIAN IN UPPER C
		 * TRIANGULAR PART AND UDIAG C EPSM --> MACHINE EPSILON C SX(N) -->
		 * DIAGONAL SCALING MATRIX FOR X C UDIAG(N) <-- ON EXIT: CONTAINS
		 * DIAGONAL OF HESSIAN C C INTERNAL VARIABLES C ------------------ C TOL
		 * TOLERANCE C DIAGMN MINIMUM ELEMENT ON DIAGONAL OF A C DIAGMX MAXIMUM
		 * ELEMENT ON DIAGONAL OF A C OFFMAX MAXIMUM OFF-DIAGONAL ELEMENT OF A C
		 * OFFROW SUM OF OFF-DIAGONAL ELEMENTS IN A ROW OF A C EVMIN MINIMUM
		 * EIGENVALUE OF A C EVMAX MAXIMUM EIGENVALUE OF A C C DESCRIPTION C
		 * ----------- C 1. IF "A" HAS ANY NEGATIVE DIAGONAL ELEMENTS, THEN
		 * CHOOSE MU>0 C SUCH THAT THE DIAGONAL OF A:=A+MU*I IS ALL POSITIVE C
		 * WITH THE RATIO OF ITS SMALLEST TO LARGEST ELEMENT ON THE C ORDER OF
		 * SQRT(EPSM). C C 2. "A" UNDERGOES A PERTURBED CHOLESKY DECOMPOSITION
		 * WHICH C RESULTS IN AN LL+ DECOMPOSITION OF A+D, WHERE D IS A C
		 * NON-NEGATIVE DIAGONAL MATRIX WHICH IS IMPLICITLY ADDED TO C "A"
		 * DURING THE DECOMPOSITION IF "A" IS NOT POSITIVE DEFINITE. C "A" IS
		 * RETAINED AND NOT CHANGED DURING THIS PROCESS BY C COPYING L INTO THE
		 * UPPER TRIANGULAR PART OF "A" AND THE C DIAGONAL INTO UDIAG. THEN THE
		 * CHOLESKY DECOMPOSITION ROUTINE C IS CALLED. ON RETURN, ADDMAX
		 * CONTAINS MAXIMUM ELEMENT OF D. C C 3. IF ADDMAX=0, "A" WAS POSITIVE
		 * DEFINITE GOING INTO STEP 2 C AND RETURN IS MADE TO CALLING PROGRAM.
		 * OTHERWISE, C THE MINIMUM NUMBER SDD WHICH MUST BE ADDED TO THE C
		 * DIAGONAL OF A TO MAKE IT SAFELY STRICTLY DIAGONALLY DOMINANT C IS
		 * CALCULATED. SINCE A+ADDMAX*I AND A+SDD*I ARE SAFELY C POSITIVE
		 * DEFINITE, CHOOSE MU=MIN(ADDMAX,SDD) AND DECOMPOSE C A+MU*I TO OBTAIN
		 * L. C
		 */

		int i, j, im1, jm1;
		double tol, diagmx, diagmn, posmax, amu, offmax, evmin, evmax, offrow, sdd;

		double addmax[] = new double[2];

		// SCALE HESSIAN
		// PRE- AND POST- MULTIPLY "A" BY INV(SX)

		for (j = 1; j <= n; j++) {

			for (i = j; i <= n; i++) {

				a[i][j] /= (sx[i] * sx[j]);

			}

		}

		// STEP1
		// -----
		// NOTE: IF A DIFFERENT TOLERANCE IS DESIRED THROUGHOUT THIS
		// ALGORITHM, CHANGE TOLERANCE HERE:

		tol = Math.sqrt(epsm);

		diagmx = a[1][1];
		diagmn = a[1][1];

		for (i = 2; i <= n; i++) {

			if (a[i][i] < diagmn)
				diagmn = a[i][i];
			if (a[i][i] > diagmx)
				diagmx = a[i][i];

		}

		posmax = Math.max(diagmx, 0.0);

		// DIAGMN .LE. 0

		if (diagmn <= posmax * tol) {

			amu = tol * (posmax - diagmn) - diagmn;

			if (amu == 0.0) {

				// FIND LARGEST OFF-DIAGONAL ELEMENT OF A

				offmax = 0.0;

				for (i = 2; i <= n; i++) {

					im1 = i - 1;

					for (j = 1; j <= im1; j++) {

						if (Math.abs(a[i][j]) > offmax)
							offmax = Math.abs(a[i][j]);

					}

				}

				amu = offmax;

				if (amu == 0.0) {

					amu = 1.0;

				} else {

					amu *= 1.0 + tol;

				}

			}

			// A = A + MU*I

			for (i = 1; i <= n; i++) {

				a[i][i] += amu;

			}

			diagmx += amu;

		}

		// STEP2
		// -----
		// COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART
		// AND DIAGONAL OF "A" TO UDIAG

		for (j = 1; j <= n; j++) {

			udiag[j] = a[j][j];

			for (i = j + 1; i <= n; i++) {

				a[j][i] = a[i][j];

			}

		}

		Uncmin_f77.choldc_f77(n, a, diagmx, tol, addmax);

		// STEP3
		// -----
		// IF ADDMAX=0, "A" WAS POSITIVE DEFINITE GOING INTO STEP 2,
		// THE LL+ DECOMPOSITION HAS BEEN DONE, AND WE RETURN.
		// OTHERWISE, ADDMAX>0. PERTURB "A" SO THAT IT IS SAFELY
		// DIAGONALLY DOMINANT AND FIND LL+ DECOMPOSITION

		if (addmax[1] > 0.0) {

			// RESTORE ORIGINAL "A" (LOWER TRIANGULAR PART AND DIAGONAL)

			for (j = 1; j <= n; j++) {

				a[j][j] = udiag[j];

				for (i = j + 1; i <= n; i++) {

					a[i][j] = a[j][i];

				}

			}

			// FIND SDD SUCH THAT A+SDD*I IS SAFELY POSITIVE DEFINITE
			// NOTE: EVMIN<0 SINCE A IS NOT POSITIVE DEFINITE;

			evmin = 0.0;
			evmax = a[1][1];

			for (i = 1; i <= n; i++) {

				offrow = 0.0;
				im1 = i - 1;

				for (j = 1; j <= im1; j++) {

					offrow += Math.abs(a[i][j]);

				}

				for (j = i + 1; j <= n; j++) {

					offrow += Math.abs(a[j][i]);

				}

				evmin = Math.min(evmin, a[i][i] - offrow);
				evmax = Math.max(evmax, a[i][i] + offrow);

			}

			sdd = tol * (evmax - evmin) - evmin;

			// PERTURB "A" AND DECOMPOSE AGAIN

			amu = Math.min(sdd, addmax[1]);

			for (i = 1; i <= n; i++) {

				a[i][i] += amu;
				udiag[i] = a[i][i];

			}

			// "A" NOW GUARANTEED SAFELY POSITIVE DEFINITE

			Uncmin_f77.choldc_f77(n, a, 0.0, tol, addmax);

		}

		// UNSCALE HESSIAN AND CHOLESKY DECOMPOSITION MATRIX

		for (j = 1; j <= n; j++) {

			for (i = j; i <= n; i++) {

				a[i][j] *= sx[i];

			}

			jm1 = j - 1;

			for (i = 1; i <= jm1; i++) {

				a[i][j] *= sx[i] * sx[j];

			}

			udiag[j] *= sx[j] * sx[j];

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The choldc_f77 method finds "THE PERTURBED L(L-TRANSPOSE) [WRITTEN LL+]
	 * DECOMPOSITION OF A+D, WHERE D IS A NON-NEGATIVE DIAGONAL MATRIX ADDED TO
	 * A IF NECESSARY TO ALLOW THE CHOLESKY DECOMPOSITION TO CONTINUE."
	 * 
	 *Translated by Steve Verrill, April 15, 1998.
	 * 
	 *@param n
	 *            Dimension of the problem
	 *@param a
	 *            On entry: matrix for which to find the perturbed Cholesky
	 *            decomposition On exit: contains L of the LL+ decomposition in
	 *            lower triangle
	 *@param diagmx
	 *            Maximum diagonal element of "A"
	 *@param tol
	 *            Tolerance
	 *@param addmax
	 *            Maximum amount implicitly added to diagonal of "A" in forming
	 *            the Cholesky decomposition of A+D
	 * 
	 */

	public static void choldc_f77(int n, double a[][], double diagmx,
			double tol, double addmax[]) {

		/*
		 * 
		 * Here is a copy of the choldc FORTRAN documentation:
		 * 
		 * SUBROUTINE CHOLDC(NR,N,A,DIAGMX,TOL,ADDMAX)
		 * 
		 * C C PURPOSE C ------- C FIND THE PERTURBED L(L-TRANSPOSE) [WRITTEN
		 * LL+] DECOMPOSITION C OF A+D, WHERE D IS A NON-NEGATIVE DIAGONAL
		 * MATRIX ADDED TO A IF C NECESSARY TO ALLOW THE CHOLESKY DECOMPOSITION
		 * TO CONTINUE. C C PARAMETERS C ---------- C NR --> ROW DIMENSION OF
		 * MATRIX C N --> DIMENSION OF PROBLEM C A(N,N) <--> ON ENTRY: MATRIX
		 * FOR WHICH TO FIND PERTURBED C CHOLESKY DECOMPOSITION C ON EXIT:
		 * CONTAINS L OF LL+ DECOMPOSITION C IN LOWER TRIANGULAR PART AND
		 * DIAGONAL OF "A" C DIAGMX --> MAXIMUM DIAGONAL ELEMENT OF "A" C TOL
		 * --> TOLERANCE C ADDMAX <-- MAXIMUM AMOUNT IMPLICITLY ADDED TO
		 * DIAGONAL OF "A" C IN FORMING THE CHOLESKY DECOMPOSITION OF A+D C
		 * INTERNAL VARIABLES C ------------------ C AMINL SMALLEST ELEMENT
		 * ALLOWED ON DIAGONAL OF L C AMNLSQ =AMINL**2 C OFFMAX MAXIMUM
		 * OFF-DIAGONAL ELEMENT IN COLUMN OF A C C C DESCRIPTION C ----------- C
		 * THE NORMAL CHOLESKY DECOMPOSITION IS PERFORMED. HOWEVER, IF AT ANY C
		 * POINT THE ALGORITHM WOULD ATTEMPT TO SET L(I,I)=SQRT(TEMP) C WITH
		 * TEMP < TOL*DIAGMX, THEN L(I,I) IS SET TO SQRT(TOL*DIAGMX) C INSTEAD.
		 * THIS IS EQUIVALENT TO ADDING TOL*DIAGMX-TEMP TO A(I,I) C C
		 */

		int i, j, jm1, jp1, k;
		double aminl, amnlsq, offmax, sum, temp;

		addmax[1] = 0.0;
		aminl = Math.sqrt(diagmx * tol);
		amnlsq = aminl * aminl;

		// FORM COLUMN J OF L

		for (j = 1; j <= n; j++) {

			// FIND DIAGONAL ELEMENTS OF L

			sum = 0.0;
			jm1 = j - 1;
			jp1 = j + 1;

			for (k = 1; k <= jm1; k++) {

				sum += a[j][k] * a[j][k];

			}

			temp = a[j][j] - sum;

			if (temp >= amnlsq) {

				a[j][j] = Math.sqrt(temp);

			} else {

				// FIND MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN

				offmax = 0.0;

				for (i = jp1; i <= n; i++) {

					if (Math.abs(a[i][j]) > offmax)
						offmax = Math.abs(a[i][j]);

				}

				if (offmax <= amnlsq)
					offmax = amnlsq;

				// ADD TO DIAGONAL ELEMENT TO ALLOW CHOLESKY DECOMPOSITION TO
				// CONTINUE

				a[j][j] = Math.sqrt(offmax);
				addmax[1] = Math.max(addmax[1], offmax - temp);

			}

			// FIND I,J ELEMENT OF LOWER TRIANGULAR MATRIX

			for (i = jp1; i <= n; i++) {

				sum = 0.0;

				for (k = 1; k <= jm1; k++) {

					sum += a[i][k] * a[j][k];

				}

				a[i][j] = (a[i][j] - sum) / a[j][j];

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The dfault_f77 method sets default values for each input variable to the
	 * minimization algorithm.
	 * 
	 *Translated by Steve Verrill, August 4, 1998.
	 * 
	 *@param n
	 *            Dimension of the problem
	 *@param x
	 *            Initial estimate of the solution (to compute max step size)
	 *@param typsiz
	 *            Typical size for each component of x
	 *@param fscale
	 *            Estimate of the scale of the minimization function
	 *@param method
	 *            Algorithm to use to solve the minimization problem
	 *@param iexp
	 *            = 0 if the minimization function is not expensive to evaluate
	 *@param msg
	 *            Message to inhibit certain automatic checks and output
	 *@param ndigit
	 *            Number of good digits in the minimization function
	 *@param itnlim
	 *            Maximum number of allowable iterations
	 *@param iagflg
	 *            = 0 if an analytic gradient is not supplied
	 *@param iahflg
	 *            = 0 if an analytic Hessian is not supplied
	 *@param dlt
	 *            Trust region radius
	 *@param gradtl
	 *            Tolerance at which the gradient is considered close enough to
	 *            zero to terminate the algorithm
	 *@param stepmx
	 *            "Value of zero to trip default maximum in optchk"
	 *@param steptl
	 *            Tolerance at which successive iterates are considered close
	 *            enough to terminate the algorithm
	 * 
	 */

	public static void dfault_f77(int n, double x[], double typsiz[],
			double fscale[], int method[], int iexp[], int msg[], int ndigit[],
			int itnlim[], int iagflg[], int iahflg[], double dlt[],
			double gradtl[], double stepmx[], double steptl[]) {

		/*
		 * 
		 * Here is a copy of the dfault FORTRAN documentation:
		 * 
		 * 
		 * SUBROUTINE DFAULT(N,X,TYPSIZ,FSCALE,METHOD,IEXP,MSG,NDIGIT, +
		 * ITNLIM,IAGFLG,IAHFLG,IPR,DLT,GRADTL,STEPMX,STEPTL)
		 * 
		 * C C PURPOSE C ------- C SET DEFAULT VALUES FOR EACH INPUT VARIABLE TO
		 * C MINIMIZATION ALGORITHM. C C PARAMETERS C ---------- C N -->
		 * DIMENSION OF PROBLEM C X(N) --> INITIAL GUESS TO SOLUTION (TO COMPUTE
		 * MAX STEP SIZE) C TYPSIZ(N) <-- TYPICAL SIZE FOR EACH COMPONENT OF X C
		 * FSCALE <-- ESTIMATE OF SCALE OF MINIMIZATION FUNCTION C METHOD <--
		 * ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM C IEXP <-- =0 IF
		 * MINIMIZATION FUNCTION NOT EXPENSIVE TO EVALUATE C MSG <-- MESSAGE TO
		 * INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT C NDIGIT <-- NUMBER OF GOOD
		 * DIGITS IN MINIMIZATION FUNCTION C ITNLIM <-- MAXIMUM NUMBER OF
		 * ALLOWABLE ITERATIONS C IAGFLG <-- =0 IF ANALYTIC GRADIENT NOT
		 * SUPPLIED C IAHFLG <-- =0 IF ANALYTIC HESSIAN NOT SUPPLIED C IPR <--
		 * DEVICE TO WHICH TO SEND OUTPUT C DLT <-- TRUST REGION RADIUS C GRADTL
		 * <-- TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE ENOUGH C TO ZERO TO
		 * TERMINATE ALGORITHM C STEPMX <-- VALUE OF ZERO TO TRIP DEFAULT
		 * MAXIMUM IN OPTCHK C STEPTL <-- TOLERANCE AT WHICH SUCCESSIVE ITERATES
		 * CONSIDERED C CLOSE ENOUGH TO TERMINATE ALGORITHM C
		 */

		int i;
		double epsm;

		// SET TYPICAL SIZE OF X AND MINIMIZATION FUNCTION

		for (i = 1; i <= n; i++) {

			typsiz[i] = 1.0;

		}

		fscale[1] = 1.0;

		// SET TOLERANCES

		dlt[1] = -1.0;

		epsm = 1.12e-16;
		gradtl[1] = Math.pow(epsm, 1.0 / 3.0);
		steptl[1] = Math.sqrt(epsm);

		stepmx[1] = 0.0;

		// SET FLAGS

		method[1] = 1;
		iexp[1] = 1;
		msg[1] = 0;
		ndigit[1] = -1;
		itnlim[1] = 150;
		iagflg[1] = 0;
		iahflg[1] = 0;

	}

	/**
	 * 
	 *<p>
	 * The dogdrv_f77 method finds the next Newton iterate (xpls) by the double
	 * dogleg method. It drives dogstp_f77.
	 * 
	 *Translated by Steve Verrill, April 15, 1998.
	 * 
	 *@param n
	 *            Dimension of the problem
	 *@param x
	 *            The old iterate
	 *@param f
	 *            Function value at the old iterate
	 *@param g
	 *            Gradient or approximation at the old iterate
	 *@param a
	 *            Cholesky decomposition of Hessian in lower triangular part and
	 *            diagonal
	 *@param p
	 *            Newton step
	 *@param xpls
	 *            The new iterate
	 *@param fpls
	 *            Function value at the new iterate
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param sx
	 *            Scaling vector for x
	 *@param stepmx
	 *            Maximum allowable step size
	 *@param steptl
	 *            Relative step size at which successive iterates are considered
	 *            close enough to terminate the algorithm
	 *@param dlt
	 *            Trust region radius (value needs to be retained between
	 *            successive calls)
	 *@param iretcd
	 *            Return code: 0 --- satisfactory xpls found 1 --- failed to
	 *            find satisfactory xpls sufficently distinct from x
	 *@param mxtake
	 *            Boolean flag indicating that a step of maximum length was used
	 *            length
	 *@param sc
	 *            Workspace (current step)
	 *@param wrk1
	 *            Workspace (and place holding argument to tregup)
	 *@param wrk2
	 *            Workspace
	 *@param wrk3
	 *            Workspace
	 * 
	 * 
	 */

	public static void dogdrv_f77(int n, double x[], double f[], double g[],
			double a[][], double p[], double xpls[], double fpls[],
			Uncmin_methods minclass, double sx[], double stepmx[],
			double steptl[], double dlt[], int iretcd[], boolean mxtake[],
			double sc[], double wrk1[], double wrk2[], double wrk3[]) {

		/*
		 * 
		 * Here is a copy of the dogdrv FORTRAN documentation:
		 * 
		 * SUBROUTINE DOGDRV(NR,N,X,F,G,A,P,XPLS,FPLS,FCN,SX,STEPMX, +
		 * STEPTL,DLT,IRETCD,MXTAKE,SC,WRK1,WRK2,WRK3,IPR)
		 * 
		 * C C PURPOSE C ------- C FIND A NEXT NEWTON ITERATE (XPLS) BY THE
		 * DOUBLE DOGLEG METHOD C C PARAMETERS C ---------- C NR --> ROW
		 * DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> OLD
		 * ITERATE X[K-1] C F --> FUNCTION VALUE AT OLD ITERATE, F(X) C G(N) -->
		 * GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE C A(N,N) --> CHOLESKY
		 * DECOMPOSITION OF HESSIAN C IN LOWER TRIANGULAR PART AND DIAGONAL C
		 * P(N) --> NEWTON STEP C XPLS(N) <-- NEW ITERATE X[K] C FPLS <--
		 * FUNCTION VALUE AT NEW ITERATE, F(XPLS) C FCN --> NAME OF SUBROUTINE
		 * TO EVALUATE FUNCTION C SX(N) --> DIAGONAL SCALING MATRIX FOR X C
		 * STEPMX --> MAXIMUM ALLOWABLE STEP SIZE C STEPTL --> RELATIVE STEP
		 * SIZE AT WHICH SUCCESSIVE ITERATES C CONSIDERED CLOSE ENOUGH TO
		 * TERMINATE ALGORITHM C DLT <--> TRUST REGION RADIUS C [RETAIN VALUE
		 * BETWEEN SUCCESSIVE CALLS] C IRETCD <-- RETURN CODE C =0 SATISFACTORY
		 * XPLS FOUND C =1 FAILED TO FIND SATISFACTORY XPLS SUFFICIENTLY C
		 * DISTINCT FROM X C MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM
		 * LENGTH USED C SC(N) --> WORKSPACE [CURRENT STEP] C WRK1(N) -->
		 * WORKSPACE (AND PLACE HOLDING ARGUMENT TO TREGUP) C WRK2(N) -->
		 * WORKSPACE C WRK3(N) --> WORKSPACE C IPR --> DEVICE TO WHICH TO SEND
		 * OUTPUT C
		 */

		int i;
		double tmp, rnwtln;
		double fplsp[] = new double[2];
		double cln[] = new double[2];
		double eta[] = new double[2];
		boolean fstdog[] = new boolean[2];
		boolean nwtake[] = new boolean[2];

		iretcd[1] = 4;
		fstdog[1] = true;
		tmp = 0.0;

		for (i = 1; i <= n; i++) {

			tmp += sx[i] * sx[i] * p[i] * p[i];

		}

		rnwtln = Math.sqrt(tmp);

		while (iretcd[1] > 1) {

			// FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM

			Uncmin_f77.dogstp_f77(n, g, a, p, sx, rnwtln, dlt, nwtake, fstdog,
					wrk1, wrk2, cln, eta, sc, stepmx);

			// CHECK NEW POINT AND UPDATE TRUST REGION

			Uncmin_f77.tregup_f77(n, x, f, g, a, minclass, sc, sx, nwtake,
					stepmx, steptl, dlt, iretcd, wrk3, fplsp, xpls, fpls,
					mxtake, 2, wrk1);

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The dogstp_f77 method finds the new step by the double dogleg appproach.
	 * 
	 *Translated by Steve Verrill, April 21, 1998.
	 * 
	 *@param n
	 *            DIMENSION OF PROBLEM
	 *@param g
	 *            GRADIENT AT CURRENT ITERATE, G(X)
	 *@param a
	 *            CHOLESKY DECOMPOSITION OF HESSIAN IN LOWER PART AND DIAGONAL
	 *@param p
	 *            NEWTON STEP
	 *@param sx
	 *            Scaling vector for x
	 *@param rnwtln
	 *            NEWTON STEP LENGTH
	 *@param dlt
	 *            TRUST REGION RADIUS
	 *@param nwtake
	 *            BOOLEAN, = true IF NEWTON STEP TAKEN
	 *@param fstdog
	 *            BOOLEAN, = true IF ON FIRST LEG OF DOGLEG
	 *@param ssd
	 *            WORKSPACE [CAUCHY STEP TO THE MINIMUM OF THE QUADRATIC MODEL
	 *            IN THE SCALED STEEPEST DESCENT DIRECTION] [RETAIN VALUE
	 *            BETWEEN SUCCESSIVE CALLS]
	 *@param v
	 *            WORKSPACE [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
	 *@param cln
	 *            CAUCHY LENGTH [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
	 *@param eta
	 *            [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
	 *@param sc
	 *            CURRENT STEP
	 *@param stepmx
	 *            MAXIMUM ALLOWABLE STEP SIZE
	 * 
	 * 
	 */

	public static void dogstp_f77(int n, double g[], double a[][], double p[],
			double sx[], double rnwtln, double dlt[], boolean nwtake[],
			boolean fstdog[], double ssd[], double v[], double cln[],
			double eta[], double sc[], double stepmx[]) {

		/*
		 * 
		 * Here is a copy of the dogstp FORTRAN documentation:
		 * 
		 * C C PURPOSE C ------- C FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM C C
		 * PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C G(N) --> GRADIENT AT CURRENT ITERATE, G(X) C
		 * A(N,N) --> CHOLESKY DECOMPOSITION OF HESSIAN IN C LOWER PART AND
		 * DIAGONAL C P(N) --> NEWTON STEP C SX(N) --> DIAGONAL SCALING MATRIX
		 * FOR X C RNWTLN --> NEWTON STEP LENGTH C DLT <--> TRUST REGION RADIUS
		 * C NWTAKE <--> BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN C FSTDOG <-->
		 * BOOLEAN, =.TRUE. IF ON FIRST LEG OF DOGLEG C SSD(N) <--> WORKSPACE
		 * [CAUCHY STEP TO THE MINIMUM OF THE C QUADRATIC MODEL IN THE SCALED
		 * STEEPEST DESCENT C DIRECTION] [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
		 * C V(N) <--> WORKSPACE [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C CLN
		 * <--> CAUCHY LENGTH C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C ETA
		 * [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C SC(N) <-- CURRENT STEP C
		 * IPR --> DEVICE TO WHICH TO SEND OUTPUT C STEPMX --> MAXIMUM ALLOWABLE
		 * STEP SIZE C C INTERNAL VARIABLES C ------------------ C CLN LENGTH OF
		 * CAUCHY STEP C
		 */

		double alpha, beta, tmp, dot1, dot2, alam;
		int i, j;

		// CAN WE TAKE NEWTON STEP

		if (rnwtln <= dlt[1]) {

			nwtake[1] = true;

			for (i = 1; i <= n; i++) {

				sc[i] = p[i];

			}

			dlt[1] = rnwtln;

		} else {

			// NEWTON STEP TOO LONG
			// CAUCHY STEP IS ON DOUBLE DOGLEG CURVE

			nwtake[1] = false;

			if (fstdog[1]) {

				// CALCULATE DOUBLE DOGLEG CURVE (SSD)

				fstdog[1] = false;
				alpha = 0.0;

				for (i = 1; i <= n; i++) {

					alpha += (g[i] * g[i]) / (sx[i] * sx[i]);

				}

				beta = 0.0;

				for (i = 1; i <= n; i++) {

					tmp = 0.0;

					for (j = i; j <= n; j++) {

						tmp += (a[j][i] * g[j]) / (sx[j] * sx[j]);

					}

					beta += tmp * tmp;

				}

				for (i = 1; i <= n; i++) {

					ssd[i] = -(alpha / beta) * g[i] / sx[i];

				}

				cln[1] = alpha * Math.sqrt(alpha) / beta;

				eta[1] = .2 + (.8 * alpha * alpha)
						/ (-beta * Blas_f77.ddot_f77(n, g, 1, p, 1));

				for (i = 1; i <= n; i++) {

					v[i] = eta[1] * sx[i] * p[i] - ssd[i];

				}

				if (dlt[1] == -1.0)
					dlt[1] = Math.min(cln[1], stepmx[1]);

			}

			if (eta[1] * rnwtln <= dlt[1]) {

				// TAKE PARTIAL STEP IN NEWTON DIRECTION

				for (i = 1; i <= n; i++) {

					sc[i] = (dlt[1] / rnwtln) * p[i];

				}

			} else {

				if (cln[1] >= dlt[1]) {

					// TAKE STEP IN STEEPEST DESCENT DIRECTION

					for (i = 1; i <= n; i++) {

						sc[i] = (dlt[1] / cln[1]) * ssd[i] / sx[i];

					}

				} else {

					// CALCULATE CONVEX COMBINATION OF SSD AND ETA*P
					// WHICH HAS SCALED LENGTH DLT

					dot1 = Blas_f77.ddot_f77(n, v, 1, ssd, 1);
					dot2 = Blas_f77.ddot_f77(n, v, 1, v, 1);

					alam = (-dot1 + Math.sqrt((dot1 * dot1) - dot2
							* (cln[1] * cln[1] - dlt[1] * dlt[1])))
							/ dot2;

					for (i = 1; i <= n; i++) {

						sc[i] = (ssd[i] + alam * v[i]) / sx[i];

					}

				}

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The forslv_f77 method solves Ax = b where A is a lower triangular matrix.
	 * 
	 *Translated by Steve Verrill, April 21, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param a
	 *            The lower triangular matrix (preserved)
	 *@param x
	 *            The solution vector
	 *@param b
	 *            The right-hand side vector
	 * 
	 * 
	 */

	public static void forslv_f77(int n, double a[][], double x[], double b[]) {

		/*
		 * 
		 * Here is a copy of the forslv FORTRAN documentation:
		 * 
		 * SUBROUTINE FORSLV(NR,N,A,X,B)
		 * 
		 * C C PURPOSE C ------- C SOLVE AX=B WHERE A IS LOWER TRIANGULAR MATRIX
		 * C C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C A(N,N) --> LOWER TRIANGULAR MATRIX (PRESERVED)
		 * C X(N) <-- SOLUTION VECTOR C B(N) --> RIGHT-HAND SIDE VECTOR C C NOTE
		 * C ---- C IF B IS NO LONGER REQUIRED BY CALLING ROUTINE, C THEN
		 * VECTORS B AND X MAY SHARE THE SAME STORAGE. C
		 */

		int i, im1, j;
		double sum;

		// SOLVE LX=B. (FORWARD SOLVE)

		x[1] = b[1] / a[1][1];

		for (i = 2; i <= n; i++) {

			sum = 0.0;
			im1 = i - 1;

			for (j = 1; j <= im1; j++) {

				sum += a[i][j] * x[j];

			}

			x[i] = (b[i] - sum) / a[i][i];

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The fstocd_f77 method finds a central difference approximation to the
	 * gradient of the function to be minimized.
	 * 
	 *Translated by Steve Verrill, April 21, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            The point at which the gradient is to be approximated
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param sx
	 *            Scaling vector for x
	 *@param rnoise
	 *            Relative noise in the function to be minimized
	 *@param g
	 *            A central difference approximation to the gradient
	 * 
	 * 
	 */

	public static void fstocd_f77(int n, double x[], Uncmin_methods minclass,
			double sx[], double rnoise, double g[]) {

		/*
		 * 
		 * Here is a copy of the fstocd FORTRAN documentation:
		 * 
		 * SUBROUTINE FSTOCD (N, X, FCN, SX, RNOISE, G)
		 * 
		 * C PURPOSE C ------- C FIND CENTRAL DIFFERENCE APPROXIMATION G TO THE
		 * FIRST DERIVATIVE C (GRADIENT) OF THE FUNCTION DEFINED BY FCN AT THE
		 * POINT X. C C PARAMETERS C ---------- C N --> DIMENSION OF PROBLEM C X
		 * --> POINT AT WHICH GRADIENT IS TO BE APPROXIMATED. C FCN --> NAME OF
		 * SUBROUTINE TO EVALUATE FUNCTION. C SX --> DIAGONAL SCALING MATRIX FOR
		 * X. C RNOISE --> RELATIVE NOISE IN FCN [F(X)]. C G <-- CENTRAL
		 * DIFFERENCE APPROXIMATION TO GRADIENT. C C
		 */

		double stepi, xtempi, fplus, fminus, xmult;
		int i;

		// FIND I-TH STEPSIZE, EVALUATE TWO NEIGHBORS IN DIRECTION OF I-TH
		// UNIT VECTOR, AND EVALUATE I-TH COMPONENT OF GRADIENT.

		xmult = Math.pow(rnoise, 1.0 / 3.0);

		for (i = 1; i <= n; i++) {

			stepi = xmult * Math.max(Math.abs(x[i]), 1.0 / sx[i]);
			xtempi = x[i];

			x[i] = xtempi + stepi;
			fplus = minclass.f_to_minimize(x);

			x[i] = xtempi - stepi;
			fminus = minclass.f_to_minimize(x);

			x[i] = xtempi;

			g[i] = (fplus - fminus) / (2.0 * stepi);

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * This version of the fstofd_f77 method finds a finite difference
	 * approximation to the Hessian.
	 * 
	 * 
	 *Translated by Steve Verrill, April 22, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param xpls
	 *            New iterate
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param fpls
	 *            fpls[1] -- fpls[n] contains the gradient of the function to
	 *            minimize
	 *@param a
	 *            "FINITE DIFFERENCE APPROXIMATION. ONLY LOWER TRIANGULAR MATRIX
	 *            AND DIAGONAL ARE RETURNED"
	 *@param sx
	 *            Scaling vector for x
	 *@param rnoise
	 *            Relative noise in the function to be minimized
	 *@param fhat
	 *            Workspace
	 * 
	 */

	public static void fstofd_f77(int n, double xpls[],
			Uncmin_methods minclass, double fpls[], double a[][], double sx[],
			double rnoise, double fhat[]) {

		/*
		 * 
		 * Here is a copy of the fstofd FORTRAN documentation. This is not
		 * entirely relevant as this Java version of fstofd only estimates the
		 * Hessian.
		 * 
		 * SUBROUTINE FSTOFD(NR,M,N,XPLS,FCN,FPLS,A,SX,RNOISE,FHAT,ICASE)
		 * 
		 * C PURPOSE C ------- C FIND FIRST ORDER FORWARD FINITE DIFFERENCE
		 * APPROXIMATION "A" TO THE C FIRST DERIVATIVE OF THE FUNCTION DEFINED
		 * BY THE SUBPROGRAM "FNAME" C EVALUATED AT THE NEW ITERATE "XPLS". C C
		 * C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE: C 1) THE FIRST
		 * DERIVATIVE (GRADIENT) OF THE OPTIMIZATION FUNCTION "FCN C ANALYTIC
		 * USER ROUTINE HAS BEEN SUPPLIED; C 2) THE SECOND DERIVATIVE (HESSIAN)
		 * OF THE OPTIMIZATION FUNCTION C IF NO ANALYTIC USER ROUTINE HAS BEEN
		 * SUPPLIED FOR THE HESSIAN BUT C ONE HAS BEEN SUPPLIED FOR THE GRADIENT
		 * ("FCN") AND IF THE C OPTIMIZATION FUNCTION IS INEXPENSIVE TO EVALUATE
		 * C C NOTE C ---- C _M=1 (OPTIMIZATION) ALGORITHM ESTIMATES THE
		 * GRADIENT OF THE FUNCTION C (FCN). FCN(X) # F: R(N)-->R(1) C _M=N
		 * (SYSTEMS) ALGORITHM ESTIMATES THE JACOBIAN OF THE FUNCTION C FCN(X) #
		 * F: R(N)-->R(N). C _M=N (OPTIMIZATION) ALGORITHM ESTIMATES THE HESSIAN
		 * OF THE OPTIMIZATIO C FUNCTION, WHERE THE HESSIAN IS THE FIRST
		 * DERIVATIVE OF "FCN" C C PARAMETERS C ---------- C NR --> ROW
		 * DIMENSION OF MATRIX C M --> NUMBER OF ROWS IN A C N --> NUMBER OF
		 * COLUMNS IN A; DIMENSION OF PROBLEM C XPLS(N) --> NEW ITERATE: X[K] C
		 * FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION C FPLS(M) --> _M=1
		 * (OPTIMIZATION) FUNCTION VALUE AT NEW ITERATE: C FCN(XPLS) C _M=N
		 * (OPTIMIZATION) VALUE OF FIRST DERIVATIVE C (GRADIENT) GIVEN BY USER
		 * FUNCTION FCN C _M=N (SYSTEMS) FUNCTION VALUE OF ASSOCIATED C
		 * MINIMIZATION FUNCTION C A(NR,N) <-- FINITE DIFFERENCE APPROXIMATION
		 * (SEE NOTE). ONLY C LOWER TRIANGULAR MATRIX AND DIAGONAL ARE RETURNED
		 * C SX(N) --> DIAGONAL SCALING MATRIX FOR X C RNOISE --> RELATIVE NOISE
		 * IN FCN [F(X)] C FHAT(M) --> WORKSPACE C ICASE --> =1 OPTIMIZATION
		 * (GRADIENT) C =2 SYSTEMS C =3 OPTIMIZATION (HESSIAN) C C INTERNAL
		 * VARIABLES C ------------------ C STEPSZ - STEPSIZE IN THE J-TH
		 * VARIABLE DIRECTION C
		 */

		double xmult, stepsz, xtmpj;
		int i, j, nm1;

		xmult = Math.sqrt(rnoise);

		for (j = 1; j <= n; j++) {

			stepsz = xmult * Math.max(Math.abs(xpls[j]), 1.0 / sx[j]);
			xtmpj = xpls[j];
			xpls[j] = xtmpj + stepsz;

			minclass.gradient(xpls, fhat);

			xpls[j] = xtmpj;

			for (i = 1; i <= n; i++) {

				a[i][j] = (fhat[i] - fpls[i]) / stepsz;

			}

		}

		nm1 = n - 1;

		for (j = 1; j <= nm1; j++) {

			for (i = j + 1; i <= n; i++) {

				a[i][j] = (a[i][j] + a[j][i]) / 2.0;

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * This version of the fstofd_f77 method finds first order finite difference
	 * approximations for gradients.
	 * 
	 *Translated by Steve Verrill, April 22, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param xpls
	 *            New iterate
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param fpls
	 *            fpls contains the value of the function to minimize at the new
	 *            iterate
	 *@param g
	 *            finite difference approximation to the gradient
	 *@param sx
	 *            Scaling vector for x
	 *@param rnoise
	 *            Relative noise in the function to be minimized
	 * 
	 */

	public static void fstofd_f77(int n, double xpls[],
			Uncmin_methods minclass, double fpls[], double g[], double sx[],
			double rnoise) {

		/*
		 * 
		 * Here is a copy of the fstofd FORTRAN documentation. It is not
		 * entirely relevant here as this particular Java method is only used
		 * for gradient calculations.
		 * 
		 * SUBROUTINE FSTOFD(NR,M,N,XPLS,FCN,FPLS,A,SX,RNOISE,FHAT,ICASE)
		 * 
		 * C PURPOSE C ------- C FIND FIRST ORDER FORWARD FINITE DIFFERENCE
		 * APPROXIMATION "A" TO THE C FIRST DERIVATIVE OF THE FUNCTION DEFINED
		 * BY THE SUBPROGRAM "FNAME" C EVALUATED AT THE NEW ITERATE "XPLS". C C
		 * C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE: C 1) THE FIRST
		 * DERIVATIVE (GRADIENT) OF THE OPTIMIZATION FUNCTION "FCN C ANALYTIC
		 * USER ROUTINE HAS BEEN SUPPLIED; C 2) THE SECOND DERIVATIVE (HESSIAN)
		 * OF THE OPTIMIZATION FUNCTION C IF NO ANALYTIC USER ROUTINE HAS BEEN
		 * SUPPLIED FOR THE HESSIAN BUT C ONE HAS BEEN SUPPLIED FOR THE GRADIENT
		 * ("FCN") AND IF THE C OPTIMIZATION FUNCTION IS INEXPENSIVE TO EVALUATE
		 * C C NOTE C ---- C _M=1 (OPTIMIZATION) ALGORITHM ESTIMATES THE
		 * GRADIENT OF THE FUNCTION C (FCN). FCN(X) # F: R(N)-->R(1) C _M=N
		 * (SYSTEMS) ALGORITHM ESTIMATES THE JACOBIAN OF THE FUNCTION C FCN(X) #
		 * F: R(N)-->R(N). C _M=N (OPTIMIZATION) ALGORITHM ESTIMATES THE HESSIAN
		 * OF THE OPTIMIZATIO C FUNCTION, WHERE THE HESSIAN IS THE FIRST
		 * DERIVATIVE OF "FCN" C C PARAMETERS C ---------- C NR --> ROW
		 * DIMENSION OF MATRIX C M --> NUMBER OF ROWS IN A C N --> NUMBER OF
		 * COLUMNS IN A; DIMENSION OF PROBLEM C XPLS(N) --> NEW ITERATE: X[K] C
		 * FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION C FPLS(M) --> _M=1
		 * (OPTIMIZATION) FUNCTION VALUE AT NEW ITERATE: C FCN(XPLS) C _M=N
		 * (OPTIMIZATION) VALUE OF FIRST DERIVATIVE C (GRADIENT) GIVEN BY USER
		 * FUNCTION FCN C _M=N (SYSTEMS) FUNCTION VALUE OF ASSOCIATED C
		 * MINIMIZATION FUNCTION C A(NR,N) <-- FINITE DIFFERENCE APPROXIMATION
		 * (SEE NOTE). ONLY C LOWER TRIANGULAR MATRIX AND DIAGONAL ARE RETURNED
		 * C SX(N) --> DIAGONAL SCALING MATRIX FOR X C RNOISE --> RELATIVE NOISE
		 * IN FCN [F(X)] C FHAT(M) --> WORKSPACE C ICASE --> =1 OPTIMIZATION
		 * (GRADIENT) C =2 SYSTEMS C =3 OPTIMIZATION (HESSIAN) C C INTERNAL
		 * VARIABLES C ------------------ C STEPSZ - STEPSIZE IN THE J-TH
		 * VARIABLE DIRECTION C
		 */

		double xmult, stepsz, xtmpj, fhat;
		int j;

		xmult = Math.sqrt(rnoise);

		// gradient

		for (j = 1; j <= n; j++) {

			stepsz = xmult * Math.max(Math.abs(xpls[j]), 1.0 / sx[j]);
			xtmpj = xpls[j];
			xpls[j] = xtmpj + stepsz;

			fhat = minclass.f_to_minimize(xpls);

			xpls[j] = xtmpj;

			g[j] = (fhat - fpls[1]) / stepsz;

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The grdchk_f77 method checks the analytic gradient supplied by the user.
	 * 
	 *Translated by Steve Verrill, April 22, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            The location at which the gradient is to be checked
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param f
	 *            Function value
	 *@param g
	 *            Analytic gradient
	 *@param typsiz
	 *            Typical size for each component of x
	 *@param sx
	 *            Scaling vector for x: sx[i] = 1.0/typsiz[i]
	 *@param fscale
	 *            Estimate of scale of f_to_minimize
	 *@param rnf
	 *            Relative noise in f_to_minimize
	 *@param analtl
	 *            Tolerance for comparison of estimated and analytical gradients
	 *@param gest
	 *            Finite difference gradient
	 * 
	 */

	public static void grdchk_f77(int n, double x[], Uncmin_methods minclass,
			double f[], double g[], double typsiz[], double sx[],
			double fscale[], double rnf, double analtl, double gest[]) {

		/*
		 * 
		 * Here is a copy of the grdchk FORTRAN documentation:
		 * 
		 * SUBROUTINE GRDCHK(N,X,FCN,F,G,TYPSIZ,SX,FSCALE,RNF, +
		 * ANALTL,WRK1,MSG,IPR)
		 * 
		 * C C PURPOSE C ------- C CHECK ANALYTIC GRADIENT AGAINST ESTIMATED
		 * GRADIENT C C PARAMETERS C ---------- C N --> DIMENSION OF PROBLEM C
		 * X(N) --> ESTIMATE TO A ROOT OF FCN C FCN --> NAME OF SUBROUTINE TO
		 * EVALUATE OPTIMIZATION FUNCTION C MUST BE DECLARED EXTERNAL IN CALLING
		 * ROUTINE C FCN: R(N) --> R(1) C F --> FUNCTION VALUE: FCN(X) C G(N)
		 * --> GRADIENT: G(X) C TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF
		 * X C SX(N) --> DIAGONAL SCALING MATRIX: SX(I)=1./TYPSIZ(I) C FSCALE
		 * --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN C RNF --> RELATIVE
		 * NOISE IN OPTIMIZATION FUNCTION FCN C ANALTL --> TOLERANCE FOR
		 * COMPARISON OF ESTIMATED AND C ANALYTICAL GRADIENTS C WRK1(N) -->
		 * WORKSPACE C MSG <-- MESSAGE OR ERROR CODE C ON OUTPUT: =-21, PROBABLE
		 * CODING ERROR OF GRADIENT C IPR --> DEVICE TO WHICH TO SEND OUTPUT C
		 */

		double gs;
		int ker, i;

		// COMPUTE FIRST ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO
		// ANALYTIC GRADIENT.

		Uncmin_f77.fstofd_f77(n, x, minclass, f, gest, sx, rnf);

		ker = 0;

		for (i = 1; i <= n; i++) {

			gs = Math.max(Math.abs(f[1]), fscale[1])
					/ Math.max(Math.abs(x[i]), typsiz[i]);

			if (Math.abs(g[i] - gest[i]) > Math.max(Math.abs(g[i]), gs)
					* analtl)
				ker = 1;

		}

		if (ker == 0)
			return;

		System.out.print("\nThere appears to be an error in the coding");
		System.out.print(" of the gradient method.\n\n\n");
		System.out.print("Component   Analytic   Finite Difference\n\n");

		for (i = 1; i <= n; i++) {

			System.out.println(i + "  " + g[i] + "  " + gest[i]);

		}

		System.exit(0);

	}

	/**
	 * 
	 *<p>
	 * The heschk_f77 method checks the analytic Hessian supplied by the user.
	 * 
	 *Translated by Steve Verrill, April 23, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            The location at which the Hessian is to be checked
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param f
	 *            Function value
	 *@param g
	 *            Gradient
	 *@param a
	 *            On exit: Hessian in lower triangle
	 *@param typsiz
	 *            Typical size for each component of x
	 *@param sx
	 *            Scaling vector for x: sx[i] = 1.0/typsiz[i]
	 *@param rnf
	 *            Relative noise in f_to_minimize
	 *@param analtl
	 *            Tolerance for comparison of estimated and analytic gradients
	 *@param iagflg
	 *            = 1 if an analytic gradient is supplied
	 *@param udiag
	 *            Workspace
	 *@param wrk1
	 *            Workspace
	 *@param wrk2
	 *            Workspace
	 */

	public static void heschk_f77(int n, double x[], Uncmin_methods minclass,
			double f[], double g[], double a[][], double typsiz[], double sx[],
			double rnf, double analtl, int iagflg[], double udiag[],
			double wrk1[], double wrk2[]) {

		/*
		 * 
		 * Here is a copy of the heschk FORTRAN documentation:
		 * 
		 * SUBROUTINE HESCHK(NR,N,X,FCN,D1FCN,D2FCN,F,G,A,TYPSIZ,SX,RNF, +
		 * ANALTL,IAGFLG,UDIAG,WRK1,WRK2,MSG,IPR)
		 * 
		 * C C PURPOSE C ------- C CHECK ANALYTIC HESSIAN AGAINST ESTIMATED
		 * HESSIAN C (THIS MAY BE DONE ONLY IF THE USER SUPPLIED ANALYTIC
		 * HESSIAN C D2FCN FILLS ONLY THE LOWER TRIANGULAR PART AND DIAGONAL OF
		 * A) C C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N
		 * --> DIMENSION OF PROBLEM C X(N) --> ESTIMATE TO A ROOT OF FCN C FCN
		 * --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION C MUST BE
		 * DECLARED EXTERNAL IN CALLING ROUTINE C FCN: R(N) --> R(1) C D1FCN -->
		 * NAME OF SUBROUTINE TO EVALUATE GRADIENT OF FCN. C MUST BE DECLARED
		 * EXTERNAL IN CALLING ROUTINE C D2FCN --> NAME OF SUBROUTINE TO
		 * EVALUATE HESSIAN OF FCN. C MUST BE DECLARED EXTERNAL IN CALLING
		 * ROUTINE C F --> FUNCTION VALUE: FCN(X) C G(N) <-- GRADIENT: G(X) C
		 * A(N,N) <-- ON EXIT: HESSIAN IN LOWER TRIANGULAR PART AND DIAG C
		 * TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF X C SX(N) -->
		 * DIAGONAL SCALING MATRIX: SX(I)=1./TYPSIZ(I) C RNF --> RELATIVE NOISE
		 * IN OPTIMIZATION FUNCTION FCN C ANALTL --> TOLERANCE FOR COMPARISON OF
		 * ESTIMATED AND C ANALYTICAL GRADIENTS C IAGFLG --> =1 IF ANALYTIC
		 * GRADIENT SUPPLIED C UDIAG(N) --> WORKSPACE C WRK1(N) --> WORKSPACE C
		 * WRK2(N) --> WORKSPACE C MSG <--> MESSAGE OR ERROR CODE C ON INPUT :
		 * IF =1XX DO NOT COMPARE ANAL + EST HESS C ON OUTPUT: =-22, PROBABLE
		 * CODING ERROR OF HESSIAN C IPR --> DEVICE TO WHICH TO SEND OUTPUT C
		 */

		int i, j, ker;
		double hs;

		// COMPUTE FINITE DIFFERENCE APPROXIMATION H TO THE HESSIAN.

		if (iagflg[1] == 1)
			Uncmin_f77.fstofd_f77(n, x, minclass, g, a, sx, rnf, wrk1);

		if (iagflg[1] != 1)
			Uncmin_f77.sndofd_f77(n, x, minclass, f, a, sx, rnf, wrk1, wrk2);

		ker = 0;

		// COPY LOWER TRIANGULAR PART OF H TO UPPER TRIANGULAR PART
		// AND DIAGONAL OF H TO UDIAG

		for (j = 1; j <= n; j++) {

			udiag[j] = a[j][j];

			for (i = j + 1; i <= n; i++) {

				a[j][i] = a[i][j];

			}

		}

		// COMPUTE ANALYTIC HESSIAN AND COMPARE TO FINITE DIFFERENCE
		// APPROXIMATION.

		minclass.hessian(x, a);

		for (j = 1; j <= n; j++) {

			hs = Math.max(Math.abs(g[j]), 1.0)
					/ Math.max(Math.abs(x[j]), typsiz[j]);

			if (Math.abs(a[j][j] - udiag[j]) > Math.max(Math.abs(udiag[j]), hs)
					* analtl)
				ker = 1;

			for (i = j + 1; i <= n; i++) {

				if (Math.abs(a[i][j] - a[j][i]) > Math.max(Math.abs(a[i][j]),
						hs)
						* analtl)
					ker = 1;

			}

		}

		if (ker == 0)
			return;

		System.out.print("\nThere appears to be an error in the coding");
		System.out.print(" of the Hessian method.\n\n\n");
		System.out.print("Row   Column   Analytic   Finite Difference\n\n");

		for (i = 1; i <= n; i++) {

			for (j = 1; j < i; j++) {

				System.out.println(i + "  " + j + "  " + a[i][j] + "  "
						+ a[j][i]);

			}

			System.out.println(i + "  " + i + "  " + a[i][i] + "  " + udiag[i]);

		}

		System.exit(0);

	}

	/**
	 * 
	 *<p>
	 * The hookdr_f77 method finds a next Newton iterate (xpls) by the
	 * More-Hebdon technique. It drives hookst_f77.
	 * 
	 *Translated by Steve Verrill, April 23, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            The old iterate
	 *@param f
	 *            The function value at the old iterate
	 *@param g
	 *            Gradient or approximation at old iterate
	 *@param a
	 *            Cholesky decomposition of Hessian in lower triangle and
	 *            diagonal. Hessian in upper triangle and udiag.
	 *@param udiag
	 *            Diagonal of Hessian in a
	 *@param p
	 *            Newton step
	 *@param xpls
	 *            New iterate
	 *@param fpls
	 *            Function value at the new iterate
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param sx
	 *            Scaling vector for x
	 *@param stepmx
	 *            Maximum allowable step size
	 *@param steptl
	 *            Relative step size at which consecutive iterates are
	 *            considered close enough to terminate the algorithm
	 *@param dlt
	 *            Trust region radius
	 *@param iretcd
	 *            Return code = 0 satisfactory xpls found = 1 failed to find
	 *            satisfactory xpls sufficiently distinct from x
	 *@param mxtake
	 *            Boolean flag indicating step of maximum length used
	 *@param amu
	 *            [Retain value between successive calls]
	 *@param dltp
	 *            [Retain value between successive calls]
	 *@param phi
	 *            [Retain value between successive calls]
	 *@param phip0
	 *            [Retain value between successive calls]
	 *@param sc
	 *            Workspace
	 *@param xplsp
	 *            Workspace
	 *@param wrk0
	 *            Workspace
	 *@param epsm
	 *            Machine epsilon
	 *@param itncnt
	 *            Iteration count
	 * 
	 */

	public static void hookdr_f77(int n, double x[], double f[], double g[],
			double a[][], double udiag[], double p[], double xpls[],
			double fpls[], Uncmin_methods minclass, double sx[],
			double stepmx[], double steptl[], double dlt[], int iretcd[],
			boolean mxtake[], double amu[], double dltp[], double phi[],
			double phip0[], double sc[], double xplsp[], double wrk0[],
			double epsm, int itncnt[]) {

		/*
		 * 
		 * Here is a copy of the hookdr FORTRAN documentation:
		 * 
		 * SUBROUTINE HOOKDR(NR,N,X,F,G,A,UDIAG,P,XPLS,FPLS,FCN,SX,STEPMX, +
		 * STEPTL,DLT,IRETCD,MXTAKE,AMU,DLTP,PHI,PHIP0, +
		 * SC,XPLSP,WRK0,EPSM,ITNCNT,IPR)
		 * 
		 * C C PURPOSE C ------- C FIND A NEXT NEWTON ITERATE (XPLS) BY THE
		 * MORE-HEBDON METHOD C C PARAMETERS C ---------- C NR --> ROW DIMENSION
		 * OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> OLD ITERATE X[K-1]
		 * C F --> FUNCTION VALUE AT OLD ITERATE, F(X) C G(N) --> GRADIENT AT
		 * OLD ITERATE, G(X), OR APPROXIMATE C A(N,N) --> CHOLESKY DECOMPOSITION
		 * OF HESSIAN IN LOWER C TRIANGULAR PART AND DIAGONAL. C HESSIAN IN
		 * UPPER TRIANGULAR PART AND UDIAG. C UDIAG(N) --> DIAGONAL OF HESSIAN
		 * IN A(.,.) C P(N) --> NEWTON STEP C XPLS(N) <-- NEW ITERATE X[K] C
		 * FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS) C FCN --> NAME OF
		 * SUBROUTINE TO EVALUATE FUNCTION C SX(N) --> DIAGONAL SCALING MATRIX
		 * FOR X C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE C STEPTL --> RELATIVE
		 * STEP SIZE AT WHICH SUCCESSIVE ITERATES C CONSIDERED CLOSE ENOUGH TO
		 * TERMINATE ALGORITHM C DLT <--> TRUST REGION RADIUS C IRETCD <--
		 * RETURN CODE C =0 SATISFACTORY XPLS FOUND C =1 FAILED TO FIND
		 * SATISFACTORY XPLS SUFFICIENTLY C DISTINCT FROM X C MXTAKE <-- BOOLEAN
		 * FLAG INDICATING STEP OF MAXIMUM LENGTH USED C AMU <--> [RETAIN VALUE
		 * BETWEEN SUCCESSIVE CALLS] C DLTP <--> [RETAIN VALUE BETWEEN
		 * SUCCESSIVE CALLS] C PHI <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
		 * C PHIP0 <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C SC(N) -->
		 * WORKSPACE C XPLSP(N) --> WORKSPACE C WRK0(N) --> WORKSPACE C EPSM -->
		 * MACHINE EPSILON C ITNCNT --> ITERATION COUNT C IPR --> DEVICE TO
		 * WHICH TO SEND OUTPUT C
		 */

		int i, j;
		boolean fstime[] = new boolean[2];
		boolean nwtake[] = new boolean[2];
		double tmp, rnwtln, alpha, beta;

		double fplsp[] = new double[2];

		iretcd[1] = 4;
		fstime[1] = true;

		tmp = 0.0;

		for (i = 1; i <= n; i++) {

			tmp += sx[i] * sx[i] * p[i] * p[i];

		}

		rnwtln = Math.sqrt(tmp);

		if (itncnt[1] == 1) {

			amu[1] = 0.0;

			// IF FIRST ITERATION AND TRUST REGION NOT PROVIDED BY USER,
			// COMPUTE INITIAL TRUST REGION.

			if (dlt[1] == -1.0) {

				alpha = 0.0;

				for (i = 1; i <= n; i++) {

					alpha += (g[i] * g[i]) / (sx[i] * sx[i]);

				}

				beta = 0.0;

				for (i = 1; i <= n; i++) {

					tmp = 0.0;

					for (j = i; j <= n; j++) {

						tmp += (a[j][i] * g[j]) / (sx[j] * sx[j]);

					}

					beta += tmp * tmp;

				}

				dlt[1] = alpha * Math.sqrt(alpha) / beta;
				dlt[1] = Math.min(dlt[1], stepmx[1]);

			}

		}

		while (iretcd[1] > 1) {

			// FIND NEW STEP BY MORE-HEBDON ALGORITHM

			Uncmin_f77.hookst_f77(n, g, a, udiag, p, sx, rnwtln, dlt, amu,
					dltp, phi, phip0, fstime, sc, nwtake, wrk0, epsm);

			dltp[1] = dlt[1];

			// CHECK NEW POINT AND UPDATE TRUST REGION

			Uncmin_f77.tregup_f77(n, x, f, g, a, minclass, sc, sx, nwtake,
					stepmx, steptl, dlt, iretcd, xplsp, fplsp, xpls, fpls,
					mxtake, 3, udiag);

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The hookst_f77 method finds a new step by the More-Hebdon algorithm. It
	 * is driven by hookdr_f77.
	 * 
	 *Translated by Steve Verrill, April 24, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param g
	 *            The gradient at the current iterate
	 *@param a
	 *            Cholesky decomposition of the Hessian in the lower triangle
	 *            and diagonal. Hessian or approximation in upper triangle (and
	 *            udiag).
	 *@udiag Diagonal of Hessian in a
	 *@param p
	 *            Newton step
	 *@param sx
	 *            Scaling vector for x
	 *@param rnwtln
	 *            Newton step length
	 *@param dlt
	 *            Trust region radius
	 *@param amu
	 *            Retain value between successive calls
	 *@param dltp
	 *            Trust region radius at last exit from this routine
	 *@param phi
	 *            Retain value between successive calls
	 *@param phip0
	 *            Retain value between successive calls
	 *@param fstime
	 *            "= true if first entry to this routine during the k-th
	 *            iteration"
	 *@param sc
	 *            Current step
	 *@param nwtake
	 *            = true if Newton step taken
	 *@param wrk0
	 *            Workspace
	 *@param epsm
	 *            Machine epsilon
	 * 
	 */

	public static void hookst_f77(int n, double g[], double a[][],
			double udiag[], double p[], double sx[], double rnwtln,
			double dlt[], double amu[], double dltp[], double phi[],
			double phip0[], boolean fstime[], double sc[], boolean nwtake[],
			double wrk0[], double epsm) {

		/*
		 * 
		 * Here is a copy of the hookst FORTRAN documentation:
		 * 
		 * SUBROUTINE HOOKST(NR,N,G,A,UDIAG,P,SX,RNWTLN,DLT,AMU, +
		 * DLTP,PHI,PHIP0,FSTIME,SC,NWTAKE,WRK0,EPSM,IPR)
		 * 
		 * C C PURPOSE C ------- C FIND NEW STEP BY MORE-HEBDON ALGORITHM C C
		 * PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C G(N) --> GRADIENT AT CURRENT ITERATE, G(X) C
		 * A(N,N) --> CHOLESKY DECOMPOSITION OF HESSIAN IN C LOWER TRIANGULAR
		 * PART AND DIAGONAL. C HESSIAN OR APPROX IN UPPER TRIANGULAR PART C
		 * UDIAG(N) --> DIAGONAL OF HESSIAN IN A(.,.) C P(N) --> NEWTON STEP C
		 * SX(N) --> DIAGONAL SCALING MATRIX FOR N C RNWTLN --> NEWTON STEP
		 * LENGTH C DLT <--> TRUST REGION RADIUS C AMU <--> [RETAIN VALUE
		 * BETWEEN SUCCESSIVE CALLS] C DLTP --> TRUST REGION RADIUS AT LAST EXIT
		 * FROM THIS ROUTINE C PHI <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
		 * C PHIP0 <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C FSTIME <-->
		 * BOOLEAN. =.TRUE. IF FIRST ENTRY TO THIS ROUTINE C DURING K-TH
		 * ITERATION C SC(N) <-- CURRENT STEP C NWTAKE <-- BOOLEAN, =.TRUE. IF
		 * NEWTON STEP TAKEN C WRK0 --> WORKSPACE C EPSM --> MACHINE EPSILON C
		 * IPR --> DEVICE TO WHICH TO SEND OUTPUT C
		 */

		double hi, alo;
		double phip, amulo, amuup, stepln;
		int i, j;
		boolean done;

		double addmax[] = new double[2];

		// HI AND ALO ARE CONSTANTS USED IN THIS ROUTINE.
		// CHANGE HERE IF OTHER VALUES ARE TO BE SUBSTITUTED.

		phip = 0.0;
		hi = 1.5;
		alo = .75;

		if (rnwtln <= hi * dlt[1]) {

			// TAKE NEWTON STEP

			nwtake[1] = true;

			for (i = 1; i <= n; i++) {

				sc[i] = p[i];

			}

			dlt[1] = Math.min(dlt[1], rnwtln);
			amu[1] = 0.0;

			return;

		} else {

			// NEWTON STEP NOT TAKEN

			nwtake[1] = false;

			if (amu[1] > 0.0) {

				amu[1] -= (phi[1] + dltp[1]) * ((dltp[1] - dlt[1]) + phi[1])
						/ (dlt[1] * phip);

			}

			phi[1] = rnwtln - dlt[1];

			if (fstime[1]) {

				for (i = 1; i <= n; i++) {

					wrk0[i] = sx[i] * sx[i] * p[i];

				}

				// SOLVE L*Y = (SX**2)*P

				Uncmin_f77.forslv_f77(n, a, wrk0, wrk0);
				phip0[1] = -Math.pow(Blas_f77.dnrm2_f77(n, wrk0, 1), 2)
						/ rnwtln;
				fstime[1] = false;

			}

			phip = phip0[1];
			amulo = -phi[1] / phip;

			amuup = 0.0;

			for (i = 1; i <= n; i++) {

				amuup += (g[i] * g[i]) / (sx[i] * sx[i]);

			}

			amuup = Math.sqrt(amuup) / dlt[1];

			done = false;

			// TEST VALUE OF AMU; GENERATE NEXT AMU IF NECESSARY

			while (!done) {

				if (amu[1] < amulo || amu[1] > amuup) {

					amu[1] = Math.max(Math.sqrt(amulo * amuup), amuup * .001);

				}

				// COPY (H,UDIAG) TO L
				// WHERE H <-- H+AMU*(SX**2) [DO NOT ACTUALLY CHANGE (H,UDIAG)]

				for (j = 1; j <= n; j++) {

					a[j][j] = udiag[j] + amu[1] * sx[j] * sx[j];

					for (i = j + 1; i <= n; i++) {

						a[i][j] = a[j][i];

					}

				}

				// FACTOR H=L(L+)

				Uncmin_f77.choldc_f77(n, a, 0.0, Math.sqrt(epsm), addmax);

				// SOLVE H*P = L(L+)*SC = -G

				for (i = 1; i <= n; i++) {

					wrk0[i] = -g[i];

				}

				Uncmin_f77.lltslv_f77(n, a, sc, wrk0);

				// RESET H. NOTE SINCE UDIAG HAS NOT BEEN DESTROYED WE NEED DO
				// NOTHING HERE. H IS IN THE UPPER PART AND IN UDIAG, STILL
				// INTACT

				stepln = 0.0;

				for (i = 1; i <= n; i++) {

					stepln += sx[i] * sx[i] * sc[i] * sc[i];

				}

				stepln = Math.sqrt(stepln);
				phi[1] = stepln - dlt[1];

				for (i = 1; i <= n; i++) {

					wrk0[i] = sx[i] * sx[i] * sc[i];

				}

				Uncmin_f77.forslv_f77(n, a, wrk0, wrk0);

				phip = -Math.pow(Blas_f77.dnrm2_f77(n, wrk0, 1), 2) / stepln;

				if ((alo * dlt[1] <= stepln && stepln <= hi * dlt[1])
						|| (amuup - amulo <= 0.0)) {

					// SC IS ACCEPTABLE HOOKSTEP

					done = true;

				} else {

					// SC NOT ACCEPTABLE HOOKSTEP. SELECT NEW AMU

					amulo = Math.max(amulo, amu[1] - (phi[1] / phip));
					if (phi[1] < 0.0)
						amuup = Math.min(amuup, amu[1]);
					amu[1] -= (stepln * phi[1]) / (dlt[1] * phip);

				}

			}

			return;

		}

	}

	/**
	 * 
	 *<p>
	 * The hsnint_f77 method provides the initial Hessian when secant updates
	 * are being used.
	 * 
	 *Translated by Steve Verrill, April 27, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param a
	 *            Initial Hessian (lower triangular matrix)
	 *@param sx
	 *            Scaling vector for x
	 *@param method
	 *            Algorithm to use to solve the minimization problem 1,2 ---
	 *            factored secant method 3 --- unfactored secant method
	 */

	public static void hsnint_f77(int n, double a[][], double sx[],
			int method[]) {

		/*
		 * 
		 * Here is a copy of the hsnint FORTRAN documentation:
		 * 
		 * SUBROUTINE HSNINT(NR,N,A,SX,METHOD)
		 * 
		 * C C PURPOSE C ------- C PROVIDE INITIAL HESSIAN WHEN USING SECANT
		 * UPDATES C C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX
		 * C N --> DIMENSION OF PROBLEM C A(N,N) <-- INITIAL HESSIAN (LOWER
		 * TRIANGULAR MATRIX) C SX(N) --> DIAGONAL SCALING MATRIX FOR X C METHOD
		 * --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM C =1,2 FACTORED
		 * SECANT METHOD USED C =3 UNFACTORED SECANT METHOD USED C
		 */

		int i, j;

		for (j = 1; j <= n; j++) {

			if (method[1] == 3) {

				a[j][j] = sx[j] * sx[j];

			} else {

				a[j][j] = sx[j];

			}

			for (i = j + 1; i <= n; i++) {

				a[i][j] = 0.0;

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The lltslv_f77 method solves Ax = b where A has the form L(L transpose)
	 * but only the lower triangular part, L, is stored.
	 * 
	 *Translated by Steve Verrill, April 27, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param a
	 *            Matrix of form L(L transpose). On return a is unchanged.
	 *@param x
	 *            The solution vector
	 *@param b
	 *            The right-hand side vector
	 * 
	 * 
	 */

	public static void lltslv_f77(int n, double a[][], double x[], double b[]) {

		/*
		 * 
		 * Here is a copy of the lltslv FORTRAN documentation:
		 * 
		 * SUBROUTINE LLTSLV(NR,N,A,X,B)
		 * 
		 * C C PURPOSE C ------- C SOLVE AX=B WHERE A HAS THE FORM
		 * L(L-TRANSPOSE) C BUT ONLY THE LOWER TRIANGULAR PART, L, IS STORED. C
		 * C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C A(N,N) --> MATRIX OF FORM L(L-TRANSPOSE). C ON
		 * RETURN A IS UNCHANGED. C X(N) <-- SOLUTION VECTOR C B(N) -->
		 * RIGHT-HAND SIDE VECTOR C C NOTE C ---- C IF B IS NOT REQUIRED BY
		 * CALLING PROGRAM, THEN C B AND X MAY SHARE THE SAME STORAGE. C
		 */

		// FORWARD SOLVE, RESULT IN X

		Uncmin_f77.forslv_f77(n, a, x, b);

		// BACK SOLVE, RESULT IN X

		Uncmin_f77.bakslv_f77(n, a, x, x);

		return;

	}

	/**
	 * 
	 *<p>
	 * The lnsrch_f77 method finds a next Newton iterate by line search.
	 * 
	 *Translated by Steve Verrill, May 15, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            Old iterate
	 *@param f
	 *            Function value at old iterate
	 *@param g
	 *            Gradient or approximation at old iterate
	 *@param p
	 *            Non-zero Newton step
	 *@param xpls
	 *            New iterate
	 *@param fpls
	 *            Function value at new iterate
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param mxtake
	 *            Boolean flag indicating whether the step of maximum length was
	 *            used
	 *@param iretcd
	 *            Return code
	 *@param stepmx
	 *            Maximum allowable step size
	 *@param steptl
	 *            Relative step size at which successive iterates are considered
	 *            close enough to terminate the algorithm
	 *@param sx
	 *            Scaling vector for x
	 * 
	 * 
	 */

	public static void lnsrch_f77(int n, double x[], double f[], double g[],
			double p[], double xpls[], double fpls[], Uncmin_methods minclass,
			boolean mxtake[], int iretcd[], double stepmx[], double steptl[],
			double sx[]) {

		/*
		 * 
		 * Here is a copy of the lnsrch FORTRAN documentation:
		 * 
		 * SUBROUTINE LNSRCH(N,X,F,G,P,XPLS,FPLS,FCN,MXTAKE, +
		 * IRETCD,STEPMX,STEPTL,SX,IPR)
		 * 
		 * C PURPOSE C ------- C FIND A NEXT NEWTON ITERATE BY LINE SEARCH. C C
		 * PARAMETERS C ---------- C N --> DIMENSION OF PROBLEM C X(N) --> OLD
		 * ITERATE: X[K-1] C F --> FUNCTION VALUE AT OLD ITERATE, F(X) C G(N)
		 * --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE C P(N) --> NON-ZERO
		 * NEWTON STEP C XPLS(N) <-- NEW ITERATE X[K] C FPLS <-- FUNCTION VALUE
		 * AT NEW ITERATE, F(XPLS) C FCN --> NAME OF SUBROUTINE TO EVALUATE
		 * FUNCTION C IRETCD <-- RETURN CODE C MXTAKE <-- BOOLEAN FLAG
		 * INDICATING STEP OF MAXIMUM LENGTH USED C STEPMX --> MAXIMUM ALLOWABLE
		 * STEP SIZE C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE
		 * ITERATES C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM C SX(N) -->
		 * DIAGONAL SCALING MATRIX FOR X C IPR --> DEVICE TO WHICH TO SEND
		 * OUTPUT C C INTERNAL VARIABLES C ------------------ C SLN NEWTON
		 * LENGTH C RLN RELATIVE LENGTH OF NEWTON STEP C
		 */

		int i;
		double tmp, sln, scl, slp, rln, rmnlmb, almbda, tlmbda;
		double t1, t2, t3, a, b, disc, pfpls, plmbda;

		pfpls = 0.0;
		plmbda = 0.0;

		mxtake[1] = false;
		iretcd[1] = 2;

		tmp = 0.0;

		for (i = 1; i <= n; i++) {

			tmp += sx[i] * sx[i] * p[i] * p[i];

		}

		sln = Math.sqrt(tmp);

		if (sln > stepmx[1]) {

			// NEWTON STEP LONGER THAN MAXIMUM ALLOWED

			scl = stepmx[1] / sln;
			Uncmin_f77.sclmul_f77(n, scl, p, p);
			sln = stepmx[1];

		}

		slp = Blas_f77.ddot_f77(n, g, 1, p, 1);
		rln = 0.0;

		for (i = 1; i <= n; i++) {

			rln = Math.max(rln, Math.abs(p[i])
					/ Math.max(Math.abs(x[i]), 1.0 / sx[i]));

		}

		rmnlmb = steptl[1] / rln;
		almbda = 1.0;

		// LOOP
		// CHECK IF NEW ITERATE SATISFACTORY. GENERATE NEW LAMBDA IF NECESSARY.

		while (iretcd[1] >= 2) {

			for (i = 1; i <= n; i++) {

				xpls[i] = x[i] + almbda * p[i];

			}

			fpls[1] = minclass.f_to_minimize(xpls);

			if (fpls[1] <= (f[1] + slp * .0001 * almbda)) {

				// SOLUTION FOUND

				iretcd[1] = 0;
				if (almbda == 1.0 && sln > .99 * stepmx[1])
					mxtake[1] = true;

			} else {

				// SOLUTION NOT (YET) FOUND

				if (almbda < rmnlmb) {

					// NO SATISFACTORY XPLS FOUND SUFFICIENTLY DISTINCT FROM X

					iretcd[1] = 1;

				} else {

					// CALCULATE NEW LAMBDA

					if (almbda == 1.0) {

						// FIRST BACKTRACK: QUADRATIC FIT

						tlmbda = -slp / (2.0 * (fpls[1] - f[1] - slp));

					} else {

						// ALL SUBSEQUENT BACKTRACKS: CUBIC FIT

						t1 = fpls[1] - f[1] - almbda * slp;
						t2 = pfpls - f[1] - plmbda * slp;
						t3 = 1.0 / (almbda - plmbda);
						a = t3
								* (t1 / (almbda * almbda) - t2
										/ (plmbda * plmbda));
						b = t3
								* (t2 * almbda / (plmbda * plmbda) - t1
										* plmbda / (almbda * almbda));
						disc = b * b - 3.0 * a * slp;

						if (disc > b * b) {

							// ONLY ONE POSITIVE CRITICAL POINT, MUST BE MINIMUM

							tlmbda = (-b + Blas_f77.sign_f77(1.0, a)
									* Math.sqrt(disc))
									/ (3.0 * a);

						} else {

							// BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM

							tlmbda = (-b - Blas_f77.sign_f77(1.0, a)
									* Math.sqrt(disc))
									/ (3.0 * a);

						}

						if (tlmbda > .5 * almbda)
							tlmbda = .5 * almbda;

					}

					plmbda = almbda;
					pfpls = fpls[1];

					if (tlmbda < almbda / 10.0) {

						almbda *= .1;

					} else {

						almbda = tlmbda;

					}

				}

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The mvmltl_f77 method computes y = Lx where L is a lower triangular
	 * matrix stored in A.
	 * 
	 *Translated by Steve Verrill, April 27, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param a
	 *            Lower triangular matrix
	 *@param x
	 *            Operand vector
	 *@param y
	 *            Result vector
	 * 
	 * 
	 */

	public static void mvmltl_f77(int n, double a[][], double x[], double y[]) {
		/*
		 * 
		 * Here is a copy of the mvmltl FORTRAN documentation:
		 * 
		 * SUBROUTINE MVMLTL(NR,N,A,X,Y) C C PURPOSE C ------- C COMPUTE Y=LX C
		 * WHERE L IS A LOWER TRIANGULAR MATRIX STORED IN A C C PARAMETERS C
		 * ---------- C NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF
		 * PROBLEM C A(N,N) --> LOWER TRIANGULAR (N*N) MATRIX C X(N) --> OPERAND
		 * VECTOR C Y(N) <-- RESULT VECTOR C C NOTE C ---- C X AND Y CANNOT
		 * SHARE STORAGE C
		 */

		double sum;
		int i, j;

		for (i = 1; i <= n; i++) {

			sum = 0.0;

			for (j = 1; j <= i; j++) {

				sum += a[i][j] * x[j];

			}

			y[i] = sum;

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The mvmlts_f77 method computes y = Ax where A is a symmetric matrix
	 * stored in its lower triangular part.
	 * 
	 *Translated by Steve Verrill, April 27, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param a
	 *            The symmetric matrix
	 *@param x
	 *            Operand vector
	 *@param y
	 *            Result vector
	 * 
	 * 
	 */

	public static void mvmlts_f77(int n, double a[][], double x[], double y[]) {
		/*
		 * 
		 * Here is a copy of the mvmlts FORTRAN documentation:
		 * 
		 * SUBROUTINE MVMLTS(NR,N,A,X,Y)
		 * 
		 * C C PURPOSE C ------- C COMPUTE Y=AX C WHERE "A" IS A SYMMETRIC (N*N)
		 * MATRIX STORED IN ITS LOWER C TRIANGULAR PART AND X,Y ARE N-VECTORS C
		 * C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C A(N,N) --> SYMMETRIC (N*N) MATRIX STORED IN C
		 * LOWER TRIANGULAR PART AND DIAGONAL C X(N) --> OPERAND VECTOR C Y(N)
		 * <-- RESULT VECTOR C C NOTE C ---- C X AND Y CANNOT SHARE STORAGE. C
		 */

		double sum;
		int i, j;

		for (i = 1; i <= n; i++) {

			sum = 0.0;

			for (j = 1; j <= i; j++) {

				sum += a[i][j] * x[j];

			}

			for (j = i + 1; j <= n; j++) {

				sum += a[j][i] * x[j];

			}

			y[i] = sum;

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The mvmltu_f77 method computes Y = (L transpose)X where L is a lower
	 * triangular matrix stored in A (L transpose is taken implicitly).
	 * 
	 *Translated by Steve Verrill, April 27, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param a
	 *            The lower triangular matrix
	 *@param x
	 *            Operand vector
	 *@param y
	 *            Result vector
	 * 
	 * 
	 */

	public static void mvmltu_f77(int n, double a[][], double x[], double y[]) {
		/*
		 * 
		 * Here is a copy of the mvmltu FORTRAN documentation:
		 * 
		 * SUBROUTINE MVMLTU(NR,N,A,X,Y)
		 * 
		 * C C PURPOSE C ------- C COMPUTE Y=(L+)X C WHERE L IS A LOWER
		 * TRIANGULAR MATRIX STORED IN A C (L-TRANSPOSE (L+) IS TAKEN
		 * IMPLICITLY) C C PARAMETERS C ---------- C NR --> ROW DIMENSION OF
		 * MATRIX C N --> DIMENSION OF PROBLEM C A(NR,1) --> LOWER TRIANGULAR
		 * (N*N) MATRIX C X(N) --> OPERAND VECTOR C Y(N) <-- RESULT VECTOR C C
		 * NOTE C ---- C X AND Y CANNOT SHARE STORAGE C
		 */

		double sum;
		int i, j;

		for (i = 1; i <= n; i++) {

			sum = 0.0;

			for (j = i; j <= n; j++) {

				sum += a[j][i] * x[j];

			}

			y[i] = sum;

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The optchk_f77 method checks the input for reasonableness.
	 * 
	 *Translated by Steve Verrill, May 12, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            On entry, estimate of the root of f_to_minimize
	 *@param typsiz
	 *            Typical size of each component of x
	 *@param sx
	 *            Scaling vector for x
	 *@param fscale
	 *            Estimate of scale of objective function
	 *@param gradtl
	 *            Tolerance at which the gradient is considered close enough to
	 *            zero to terminate the algorithm
	 *@param itnlim
	 *            Maximum number of allowable iterations
	 *@param ndigit
	 *            Number of good digits in the optimization function
	 *@param epsm
	 *            Machine epsilon
	 *@param dlt
	 *            Trust region radius
	 *@param method
	 *            Algorithm indicator
	 *@param iexp
	 *            Expense flag
	 *@param iagflg
	 *            = 1 if an analytic gradient is supplied
	 *@param iahflg
	 *            = 1 if an analytic Hessian is supplied
	 *@param stepmx
	 *            Maximum step size
	 *@param msg
	 *            Message and error code
	 * 
	 */

	public static void optchk_f77(int n, double x[], double typsiz[],
			double sx[], double fscale[], double gradtl[], int itnlim[],
			int ndigit[], double epsm, double dlt[], int method[], int iexp[],
			int iagflg[], int iahflg[], double stepmx[], int msg[]) {

		/*
		 * 
		 * Here is a copy of the optchk FORTRAN documentation:
		 * 
		 * 
		 * SUBROUTINE OPTCHK(N,X,TYPSIZ,SX,FSCALE,GRADTL,ITNLIM,NDIGIT,EPSM, +
		 * DLT,METHOD,IEXP,IAGFLG,IAHFLG,STEPMX,MSG,IPR) C C PURPOSE C ------- C
		 * CHECK INPUT FOR REASONABLENESS C C PARAMETERS C ---------- C N -->
		 * DIMENSION OF PROBLEM C X(N) --> ON ENTRY, ESTIMATE TO ROOT OF FCN C
		 * TYPSIZ(N) <--> TYPICAL SIZE OF EACH COMPONENT OF X C SX(N) <--
		 * DIAGONAL SCALING MATRIX FOR X C FSCALE <--> ESTIMATE OF SCALE OF
		 * OBJECTIVE FUNCTION FCN C GRADTL --> TOLERANCE AT WHICH GRADIENT
		 * CONSIDERED CLOSE C ENOUGH TO ZERO TO TERMINATE ALGORITHM C ITNLIM
		 * <--> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS C NDIGIT <--> NUMBER OF
		 * GOOD DIGITS IN OPTIMIZATION FUNCTION FCN C EPSM --> MACHINE EPSILON C
		 * DLT <--> TRUST REGION RADIUS C METHOD <--> ALGORITHM INDICATOR C IEXP
		 * <--> EXPENSE FLAG C IAGFLG <--> =1 IF ANALYTIC GRADIENT SUPPLIED C
		 * IAHFLG <--> =1 IF ANALYTIC HESSIAN SUPPLIED C STEPMX <--> MAXIMUM
		 * STEP SIZE C MSG <--> MESSAGE AND ERROR CODE C IPR --> DEVICE TO WHICH
		 * TO SEND OUTPUT C
		 */

		int i;
		double stpsiz;

		// CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES.
		// IF NOT, SET THEM TO DEFAULT VALUES.

		if (method[1] < 1 || method[1] > 3)
			method[1] = 1;
		if (iagflg[1] != 1)
			iagflg[1] = 0;
		if (iahflg[1] != 1)
			iahflg[1] = 0;
		if (iexp[1] != 0)
			iexp[1] = 1;

		if ((msg[1] / 2) % 2 == 1 && iagflg[1] == 0) {

			System.out
					.print("\n\nOPTCHK   User requests that analytic gradient");
			System.out.print(" be accepted as properly coded,\n");
			System.out.print("OPTCHK   msg = " + msg + ",\n");
			System.out.print("OPTCHK   but an analytic gradient is not");
			System.out.print(" supplied,\n");
			System.out.print("OPTCHK   iagflg = " + iagflg[1] + ".\n\n");

			System.exit(0);

		}

		if ((msg[1] / 4) % 2 == 1 && iahflg[1] == 0) {

			System.out
					.print("\n\nOPTCHK   User requests that analytic Hessian");
			System.out.print(" be accepted as properly coded,\n");
			System.out.print("OPTCHK   msg = " + msg + ",\n");
			System.out.print("OPTCHK   but an analytic Hessian is not");
			System.out.print(" supplied,\n");
			System.out.print("OPTCHK   iahflg = " + iahflg[1] + ".\n\n");

			System.exit(0);

		}

		// CHECK DIMENSION OF PROBLEM

		if (n <= 0) {

			System.out.print("\n\nOPTCHK   Illegal dimension, n = " + n
					+ "\n\n");

			System.exit(0);

		}

		if (n == 1 && msg[1] % 2 == 0) {

			System.out.print("\n\nOPTCHK   !!!WARNING!!!  This class is ");
			System.out.print("inefficient for problems of size 1.\n");
			System.out
					.print("OPTCHK   You might want to try Fmin instead.\n\n");

			msg[1] = -2;

		}

		// COMPUTE SCALE MATRIX

		for (i = 1; i <= n; i++) {

			if (typsiz[i] == 0)
				typsiz[i] = 1.0;
			if (typsiz[i] < 0.0)
				typsiz[i] = -typsiz[i];
			sx[i] = 1.0 / typsiz[i];

		}

		// CHECK MAXIMUM STEP SIZE

		if (stepmx[1] <= 0.0) {

			stpsiz = 0.0;

			for (i = 1; i <= n; i++) {

				stpsiz += x[i] * x[i] * sx[i] * sx[i];

			}

			stpsiz = Math.sqrt(stpsiz);

			stepmx[1] = Math.max(1000.0 * stpsiz, 1000.0);

		}

		// CHECK FUNCTION SCALE

		if (fscale[1] == 0)
			fscale[1] = 1.0;
		if (fscale[1] < 0.0)
			fscale[1] = -fscale[1];

		// CHECK GRADIENT TOLERANCE

		if (gradtl[1] < 0.0) {

			System.out.print("\n\nOPTCHK   Illegal tolerance, gradtl = "
					+ gradtl[1] + "\n\n");

			System.exit(0);

		}

		// CHECK ITERATION LIMIT

		if (itnlim[1] < 0) {

			System.out.print("\n\nOPTCHK   Illegal iteration limit,");
			System.out.print(" itnlim = " + itnlim[1] + "\n\n");

			System.exit(0);

		}

		// CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN

		if (ndigit[1] == 0) {

			System.out.print("\n\nOPTCHK   Minimization function has no good");
			System.out.print(" digits, ndigit = " + ndigit + "\n\n");

			System.exit(0);

		}

		if (ndigit[1] < 0)
			ndigit[1] = (int) (-Math.log(epsm) / Math.log(10.0));

		// CHECK TRUST REGION RADIUS

		if (dlt[1] <= 0.0)
			dlt[1] = -1.0;
		if (dlt[1] > stepmx[1])
			dlt[1] = stepmx[1];

		return;

	}

	/**
	 * 
	 *<p>
	 * The optdrv_f77 method is the driver for the nonlinear optimization
	 * problem.
	 * 
	 *Translated by Steve Verrill, May 18, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            On entry, estimate of the location of a minimum of
	 *            f_to_minimize
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param typsiz
	 *            Typical size of each component of x
	 *@param fscale
	 *            Estimate of scale of objective function
	 *@param method
	 *            Algorithm indicator 1 -- line search 2 -- double dogleg 3 --
	 *            More-Hebdon
	 *@param iexp
	 *            Expense flag. 1 -- optimization function, f_to_minimize, is
	 *            expensive to evaluate 0 -- otherwise If iexp = 1, the Hessian
	 *            will be evaluated by secant update rather than analytically or
	 *            by finite differences.
	 *@param msg
	 *            On input: (> 0) message to inhibit certain automatic checks On
	 *            output: (< 0) error code (= 0, no error)
	 *@param ndigit
	 *            Number of good digits in the optimization function
	 *@param itnlim
	 *            Maximum number of allowable iterations
	 *@param iagflg
	 *            = 1 if an analytic gradient is supplied
	 *@param iahflg
	 *            = 1 if an analytic Hessian is supplied
	 *@param dlt
	 *            Trust region radius
	 *@param gradtl
	 *            Tolerance at which the gradient is considered close enough to
	 *            zero to terminate the algorithm
	 *@param stepmx
	 *            Maximum step size
	 *@param steptl
	 *            Relative step size at which successive iterates are considered
	 *            close enough to terminate the algorithm
	 *@param xpls
	 *            On exit: xpls is a local minimum
	 *@param fpls
	 *            On exit: function value at xpls
	 *@param gpls
	 *            On exit: gradient at xpls
	 *@param itrmcd
	 *            Termination code
	 *@param a
	 *            workspace for Hessian (or its approximation) and its Cholesky
	 *            decomposition
	 *@param udiag
	 *            workspace (for diagonal of Hessian)
	 *@param g
	 *            workspace (for gradient at current iterate)
	 *@param p
	 *            workspace for step
	 *@param sx
	 *            workspace (for scaling vector)
	 *@param wrk0
	 *            workspace
	 *@param wrk1
	 *            workspace
	 *@param wrk2
	 *            workspace
	 *@param wrk3
	 *            workspace
	 * 
	 */

	public static void optdrv_f77(int n, double x[], Uncmin_methods minclass,
			double typsiz[], double fscale[], int method[], int iexp[],
			int msg[], int ndigit[], int itnlim[], int iagflg[], int iahflg[],
			double dlt[], double gradtl[], double stepmx[], double steptl[],
			double xpls[], double fpls[], double gpls[], int itrmcd[],
			double a[][], double udiag[], double g[], double p[], double sx[],
			double wrk0[], double wrk1[], double wrk2[], double wrk3[]) {

		/*
		 * 
		 * Here is a copy of the optdrv FORTRAN documentation:
		 * 
		 * 
		 * SUBROUTINE OPTDRV(NR,N,X,FCN,D1FCN,D2FCN,TYPSIZ,FSCALE, +
		 * METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR, +
		 * DLT,GRADTL,STEPMX,STEPTL, + XPLS,FPLS,GPLS,ITRMCD, +
		 * A,UDIAG,G,P,SX,WRK0,WRK1,WRK2,WRK3)
		 * 
		 * C C PURPOSE C ------- C DRIVER FOR NON-LINEAR OPTIMIZATION PROBLEM C
		 * C PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF PROBLEM C X(N) --> ON ENTRY: ESTIMATE TO A ROOT OF FCN C
		 * FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION C MUST
		 * BE DECLARED EXTERNAL IN CALLING ROUTINE C FCN: R(N) --> R(1) C D1FCN
		 * --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE GRADIENT C OF FCN. MUST
		 * BE DECLARED EXTERNAL IN CALLING ROUTINE C D2FCN --> (OPTIONAL) NAME
		 * OF SUBROUTINE TO EVALUATE HESSIAN OF C OF FCN. MUST BE DECLARED
		 * EXTERNAL IN CALLING ROUTINE C TYPSIZ(N) --> TYPICAL SIZE FOR EACH
		 * COMPONENT OF X C FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION C
		 * METHOD --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM C =1 LINE
		 * SEARCH C =2 DOUBLE DOGLEG C =3 MORE-HEBDON C IEXP --> =1 IF
		 * OPTIMIZATION FUNCTION FCN IS EXPENSIVE TO C EVALUATE, =0 OTHERWISE.
		 * IF SET THEN HESSIAN WILL C BE EVALUATED BY SECANT UPDATE INSTEAD OF C
		 * ANALYTICALLY OR BY FINITE DIFFERENCES C MSG <--> ON INPUT: (.GT.0)
		 * MESSAGE TO INHIBIT CERTAIN C AUTOMATIC CHECKS C ON OUTPUT: (.LT.0)
		 * ERROR CODE; =0 NO ERROR C NDIGIT --> NUMBER OF GOOD DIGITS IN
		 * OPTIMIZATION FUNCTION FCN C ITNLIM --> MAXIMUM NUMBER OF ALLOWABLE
		 * ITERATIONS C IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED C IAHFLG -->
		 * =1 IF ANALYTIC HESSIAN SUPPLIED C IPR --> DEVICE TO WHICH TO SEND
		 * OUTPUT C DLT --> TRUST REGION RADIUS C GRADTL --> TOLERANCE AT WHICH
		 * GRADIENT CONSIDERED CLOSE C ENOUGH TO ZERO TO TERMINATE ALGORITHM C
		 * STEPMX --> MAXIMUM ALLOWABLE STEP SIZE C STEPTL --> RELATIVE STEP
		 * SIZE AT WHICH SUCCESSIVE ITERATES C CONSIDERED CLOSE ENOUGH TO
		 * TERMINATE ALGORITHM C XPLS(N) <--> ON EXIT: XPLS IS LOCAL MINIMUM C
		 * FPLS <--> ON EXIT: FUNCTION VALUE AT SOLUTION, XPLS C GPLS(N) <--> ON
		 * EXIT: GRADIENT AT SOLUTION XPLS C ITRMCD <-- TERMINATION CODE C
		 * A(N,N) --> WORKSPACE FOR HESSIAN (OR ESTIMATE) C AND ITS CHOLESKY
		 * DECOMPOSITION C UDIAG(N) --> WORKSPACE [FOR DIAGONAL OF HESSIAN] C
		 * G(N) --> WORKSPACE (FOR GRADIENT AT CURRENT ITERATE) C P(N) -->
		 * WORKSPACE FOR STEP C SX(N) --> WORKSPACE (FOR DIAGONAL SCALING
		 * MATRIX) C WRK0(N) --> WORKSPACE C WRK1(N) --> WORKSPACE C WRK2(N) -->
		 * WORKSPACE C WRK3(N) --> WORKSPACE C C C INTERNAL VARIABLES C
		 * ------------------ C ANALTL TOLERANCE FOR COMPARISON OF ESTIMATED AND
		 * C ANALYTICAL GRADIENTS AND HESSIANS C EPSM MACHINE EPSILON C F
		 * FUNCTION VALUE: FCN(X) C ITNCNT CURRENT ITERATION, K C RNF RELATIVE
		 * NOISE IN OPTIMIZATION FUNCTION FCN. C NOISE=10.**(-NDIGIT) C
		 */

		boolean noupdt[] = new boolean[2];
		boolean mxtake[] = new boolean[2];

		int i, j, num5, remain, ilow, ihigh;
		int icscmx[] = new int[2];
		int iretcd[] = new int[2];
		int itncnt[] = new int[2];

		double rnf, analtl, dltsav, amusav, dlpsav, phisav, phpsav;

		double f[] = new double[2];

		double amu[] = new double[2];
		double dltp[] = new double[2];
		double phi[] = new double[2];
		double phip0[] = new double[2];

		double epsm;

		dltsav = amusav = dlpsav = phisav = phpsav = 0.0;

		// INITIALIZATION

		for (i = 1; i <= n; i++) {

			p[i] = 0.0;

		}

		itncnt[1] = 0;
		iretcd[1] = -1;

		epsm = 1.12e-16;

		Uncmin_f77.optchk_f77(n, x, typsiz, sx, fscale, gradtl, itnlim, ndigit,
				epsm, dlt, method, iexp, iagflg, iahflg, stepmx, msg);

		// removed because don't want to stop if msg = -2
		// and other negative messages are now handled where they are generated
		// if (msg[1] < 0) return;

		rnf = Math.max(Math.pow(10.0, -ndigit[1]), epsm);
		analtl = Math.max(.01, Math.sqrt(rnf));

		if ((msg[1] / 8) % 2 != 1) {

			num5 = n / 5;
			remain = n % 5;

			System.out.print("\n\nOPTDRV          Typical x\n\n");

			ilow = -4;
			ihigh = 0;

			for (i = 1; i <= num5; i++) {

				ilow += 5;
				ihigh += 5;

				System.out.print(ilow + "--" + ihigh + "     ");

				for (j = 1; j <= 5; j++) {

					System.out.print(typsiz[ilow + j - 1] + "  ");

				}

				System.out.print("\n");

			}

			ilow += 5;
			ihigh = ilow + remain - 1;

			System.out.print(ilow + "--" + ihigh + "     ");

			for (j = 1; j <= remain; j++) {

				System.out.print(typsiz[ilow + j - 1] + "  ");

			}

			System.out.print("\n");

			System.out.print("\n\nOPTDRV      Scaling vector for x\n\n");

			ilow = -4;
			ihigh = 0;

			for (i = 1; i <= num5; i++) {

				ilow += 5;
				ihigh += 5;

				System.out.print(ilow + "--" + ihigh + "     ");

				for (j = 1; j <= 5; j++) {

					System.out.print(sx[ilow + j - 1] + "  ");

				}

				System.out.print("\n");

			}

			ilow += 5;
			ihigh = ilow + remain - 1;

			System.out.print(ilow + "--" + ihigh + "     ");

			for (j = 1; j <= remain; j++) {

				System.out.print(sx[ilow + j - 1] + "  ");

			}

			System.out.print("\n");

			System.out.println("\n\nOPTDRV      Typical f = " + fscale[1]);
			System.out.print("OPTDRV      Number of good digits in");
			System.out.println(" f_to_minimize = " + ndigit[1]);
			System.out.print("OPTDRV      Gradient flag");
			System.out.println(" = " + iagflg[1]);
			System.out.print("OPTDRV      Hessian flag");
			System.out.println(" = " + iahflg[1]);
			System.out.print("OPTDRV      Expensive function calculation flag");
			System.out.println(" = " + iexp[1]);
			System.out.print("OPTDRV      Method to use");
			System.out.println(" = " + method[1]);
			System.out.print("OPTDRV      Iteration limit");
			System.out.println(" = " + itnlim[1]);
			System.out.print("OPTDRV      Machine epsilon");
			System.out.println(" = " + epsm);
			System.out.print("OPTDRV      Maximum step size");
			System.out.println(" = " + stepmx[1]);
			System.out.print("OPTDRV      Step tolerance");
			System.out.println(" = " + steptl[1]);
			System.out.print("OPTDRV      Gradient tolerance");
			System.out.println(" = " + gradtl[1]);
			System.out.print("OPTDRV      Trust region radius");
			System.out.println(" = " + dlt[1]);
			System.out.print("OPTDRV      Relative noise in");
			System.out.println(" f_to_minimize = " + rnf);
			System.out.print("OPTDRV      Analytical fd tolerance");
			System.out.println(" = " + analtl);

		}

		// EVALUATE FCN(X)

		f[1] = minclass.f_to_minimize(x);

		// EVALUATE ANALYTIC OR FINITE DIFFERENCE GRADIENT AND CHECK ANALYTIC
		// GRADIENT, IF REQUESTED.

		if (iagflg[1] == 0) {

			Uncmin_f77.fstofd_f77(n, x, minclass, f, g, sx, rnf);

		} else {

			minclass.gradient(x, g);

			if ((msg[1] / 2) % 2 == 0) {

				Uncmin_f77.grdchk_f77(n, x, minclass, f, g, typsiz, sx, fscale,
						rnf, analtl, wrk1);

			}

		}

		Uncmin_f77.optstp_f77(n, x, f, g, wrk1, itncnt, icscmx, itrmcd, gradtl,
				steptl, sx, fscale, itnlim, iretcd, mxtake, msg);

		if (itrmcd[1] != 0) {

			fpls[1] = f[1];

			for (i = 1; i <= n; i++) {

				xpls[i] = x[i];
				gpls[i] = g[i];

			}

		} else {

			if (iexp[1] == 1) {

				// IF OPTIMIZATION FUNCTION EXPENSIVE TO EVALUATE (IEXP=1), THEN
				// HESSIAN WILL BE OBTAINED BY SECANT UPDATES. GET INITIAL
				// HESSIAN.

				Uncmin_f77.hsnint_f77(n, a, sx, method);

			} else {

				// EVALUATE ANALYTIC OR FINITE DIFFERENCE HESSIAN AND CHECK
				// ANALYTIC
				// HESSIAN IF REQUESTED (ONLY IF USER-SUPPLIED ANALYTIC HESSIAN
				// ROUTINE minclass.hessian FILLS ONLY LOWER TRIANGULAR PART AND
				// DIAGONAL OF A).

				if (iahflg[1] == 0) {

					if (iagflg[1] == 1) {

						Uncmin_f77.fstofd_f77(n, x, minclass, g, a, sx, rnf,
								wrk1);

					} else {

						Uncmin_f77.sndofd_f77(n, x, minclass, f, a, sx, rnf,
								wrk1, wrk2);

					}

				} else {

					if ((msg[1] / 4) % 2 == 0) {// Initially the condition was:
												// (msg[1]/4)%2 == 1

						minclass.hessian(x, a);

					} else {
						Uncmin_f77.heschk_f77(n, x, minclass, f, g, a, typsiz,
								sx, rnf, analtl, iagflg, udiag, wrk1, wrk2);

						// HESCHK EVALUATES minclass.hessian AND CHECKS IT
						// AGAINST THE FINITE
						// DIFFERENCE HESSIAN WHICH IT CALCULATES BY CALLING
						// FSTOFD
						// (IF IAGFLG .EQ. 1) OR SNDOFD (OTHERWISE).

					}

				}

			}

			if ((msg[1] / 8) % 2 == 0)
				Uncmin_f77.result_f77(n, x, f, g, a, p, itncnt, 1);

			// ITERATION

			while (itrmcd[1] == 0) {

				itncnt[1]++;

				// FIND PERTURBED LOCAL MODEL HESSIAN AND ITS LL+ DECOMPOSITION
				// (SKIP THIS STEP IF LINE SEARCH OR DOGSTEP TECHNIQUES BEING
				// USED WITH
				// SECANT UPDATES. CHOLESKY DECOMPOSITION L ALREADY OBTAINED
				// FROM
				// SECFAC.)

				if (iexp[1] != 1 || method[1] == 3) {

					Uncmin_f77.chlhsn_f77(n, a, epsm, sx, udiag);

				}

				// SOLVE FOR NEWTON STEP: AP=-G

				for (i = 1; i <= n; i++) {

					wrk1[i] = -g[i];

				}

				Uncmin_f77.lltslv_f77(n, a, p, wrk1);

				// DECIDE WHETHER TO ACCEPT NEWTON STEP XPLS=X + P
				// OR TO CHOOSE XPLS BY A GLOBAL STRATEGY.

				if (iagflg[1] == 0 && method[1] != 1) {

					dltsav = dlt[1];

					if (method[1] != 2) {

						amusav = amu[1];
						dlpsav = dltp[1];
						phisav = phi[1];
						phpsav = phip0[1];

					}

				}

				if (method[1] == 1) {

					Uncmin_f77.lnsrch_f77(n, x, f, g, p, xpls, fpls, minclass,
							mxtake, iretcd, stepmx, steptl, sx);

				} else if (method[1] == 2) {

					Uncmin_f77.dogdrv_f77(n, x, f, g, a, p, xpls, fpls,
							minclass, sx, stepmx, steptl, dlt, iretcd, mxtake,
							wrk0, wrk1, wrk2, wrk3);

				} else {

					Uncmin_f77.hookdr_f77(n, x, f, g, a, udiag, p, xpls, fpls,
							minclass, sx, stepmx, steptl, dlt, iretcd, mxtake,
							amu, dltp, phi, phip0, wrk0, wrk1, wrk2, epsm,
							itncnt);

				}

				// IF COULD NOT FIND SATISFACTORY STEP AND FORWARD DIFFERENCE
				// GRADIENT WAS USED, RETRY USING CENTRAL DIFFERENCE GRADIENT.

				if (iretcd[1] == 1 && iagflg[1] == 0) {

					// SET IAGFLG FOR CENTRAL DIFFERENCES

					iagflg[1] = -1;

					System.out
							.print("\nOPTDRV      Shift from forward to central");
					System.out.print(" differences");
					System.out.print("\nOPTDRV      in iteration " + itncnt[1]);
					System.out.print("\n");

					Uncmin_f77.fstocd_f77(n, x, minclass, sx, rnf, g);

					if (method[1] == 1) {

						// SOLVE FOR NEWTON STEP: AP=-G

						for (i = 1; i <= n; i++) {

							wrk1[i] = -g[i];

						}

						Uncmin_f77.lltslv_f77(n, a, p, wrk1);

						Uncmin_f77.lnsrch_f77(n, x, f, g, p, xpls, fpls,
								minclass, mxtake, iretcd, stepmx, steptl, sx);

					} else {

						dlt[1] = dltsav;
						if (method[1] == 2) {

							// SOLVE FOR NEWTON STEP: AP=-G

							for (i = 1; i <= n; i++) {

								wrk1[i] = -g[i];

							}

							Uncmin_f77.lltslv_f77(n, a, p, wrk1);

							Uncmin_f77.dogdrv_f77(n, x, f, g, a, p, xpls, fpls,
									minclass, sx, stepmx, steptl, dlt, iretcd,
									mxtake, wrk0, wrk1, wrk2, wrk3);

						} else {

							amu[1] = amusav;
							dltp[1] = dlpsav;
							phi[1] = phisav;
							phip0[1] = phpsav;

							Uncmin_f77.chlhsn_f77(n, a, epsm, sx, udiag);

							// SOLVE FOR NEWTON STEP: AP=-G

							for (i = 1; i <= n; i++) {

								wrk1[i] = -g[i];

							}

							Uncmin_f77.lltslv_f77(n, a, p, wrk1);

							Uncmin_f77.hookdr_f77(n, x, f, g, a, udiag, p,
									xpls, fpls, minclass, sx, stepmx, steptl,
									dlt, iretcd, mxtake, amu, dltp, phi, phip0,
									wrk0, wrk1, wrk2, epsm, itncnt);

						}

					}

				}

				// CALCULATE STEP FOR OUTPUT

				for (i = 1; i <= n; i++) {

					p[i] = xpls[i] - x[i];

				}

				// CALCULATE GRADIENT AT XPLS

				if (iagflg[1] == -1) {

					// CENTRAL DIFFERENCE GRADIENT

					Uncmin_f77.fstocd_f77(n, xpls, minclass, sx, rnf, gpls);

				} else if (iagflg[1] == 0) {

					// FORWARD DIFFERENCE GRADIENT

					Uncmin_f77.fstofd_f77(n, xpls, minclass, fpls, gpls, sx,
							rnf);

				} else {

					// ANALYTIC GRADIENT

					minclass.gradient(xpls, gpls);

				}

				// CHECK WHETHER STOPPING CRITERIA SATISFIED

				Uncmin_f77.optstp_f77(n, xpls, fpls, gpls, x, itncnt, icscmx,
						itrmcd, gradtl, steptl, sx, fscale, itnlim, iretcd,
						mxtake, msg);

				if (itrmcd[1] == 0) {

					// EVALUATE HESSIAN AT XPLS

					if (iexp[1] != 0) {

						if (method[1] == 3) {

							Uncmin_f77.secunf_f77(n, x, g, a, udiag, xpls,
									gpls, epsm, itncnt, rnf, iagflg, noupdt,
									wrk1, wrk2, wrk3);

						} else {

							Uncmin_f77.secfac_f77(n, x, g, a, xpls, gpls, epsm,
									itncnt, rnf, iagflg, noupdt, wrk0, wrk1,
									wrk2, wrk3);

						}

					} else {

						if (iahflg[1] == 1) {

							minclass.hessian(xpls, a);

						} else {

							if (iagflg[1] == 1) {

								Uncmin_f77.fstofd_f77(n, xpls, minclass, gpls,
										a, sx, rnf, wrk1);

							} else {

								Uncmin_f77.sndofd_f77(n, xpls, minclass, fpls,
										a, sx, rnf, wrk1, wrk2);

							}

						}

					}

					if ((msg[1] / 16) % 2 == 1) {

						Uncmin_f77.result_f77(n, xpls, fpls, gpls, a, p,
								itncnt, 1);

					}

					// X <-- XPLS AND G <-- GPLS AND F <-- FPLS

					f[1] = fpls[1];

					for (i = 1; i <= n; i++) {

						x[i] = xpls[i];
						g[i] = gpls[i];

					}

				}

			}

			// TERMINATION
			// -----------
			// RESET XPLS,FPLS,GPLS, IF PREVIOUS ITERATE SOLUTION
			//

			if (itrmcd[1] == 3) {

				fpls[1] = f[1];

				for (i = 1; i <= n; i++) {

					xpls[i] = x[i];
					gpls[i] = g[i];

				}

			}

		}

		// PRINT RESULTS

		if ((msg[1] / 8) % 2 == 0) {

			Uncmin_f77.result_f77(n, xpls, fpls, gpls, a, p, itncnt, 0);

		}

		msg[1] = 0;
		return;

	}

	/**
	 * 
	 *<p>
	 * The optstp_f77 method determines whether the algorithm should terminate
	 * due to any of the following: 1) problem solved within user tolerance 2)
	 * convergence within user tolerance 3) iteration limit reached 4)
	 * divergence or too restrictive maximum step (stepmx) suspected
	 * 
	 *Translated by Steve Verrill, May 12, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param xpls
	 *            New iterate
	 *@param fpls
	 *            Function value at new iterate
	 *@param gpls
	 *            Gradient or approximation at new iterate
	 *@param x
	 *            Old iterate
	 *@param itncnt
	 *            Current iteration
	 *@param icscmx
	 *            Number of consecutive steps >= stepmx (retain between
	 *            successive calls)
	 *@param itrmcd
	 *            Termination code
	 *@param gradtl
	 *            Tolerance at which the relative gradient is considered close
	 *            enough to zero to terminate the algorithm
	 *@param steptl
	 *            Relative step size at which successive iterates are considered
	 *            close enough to terminate the algorithm
	 *@param sx
	 *            Scaling vector for x
	 *@param fscale
	 *            Estimate of the scale of the objective function
	 *@param itnlim
	 *            Maximum number of allowable iterations
	 *@param iretcd
	 *            Return code
	 *@param mxtake
	 *            Boolean flag indicating step of maximum length was used
	 *@param msg
	 *            If msg includes a term 8, suppress output
	 * 
	 */

	public static void optstp_f77(int n, double xpls[], double fpls[],
			double gpls[], double x[], int itncnt[], int icscmx[],
			int itrmcd[], double gradtl[], double steptl[], double sx[],
			double fscale[], int itnlim[], int iretcd[], boolean mxtake[],
			int msg[]) {
		/*
		 * 
		 * Here is a copy of the optstp FORTRAN documentation:
		 * 
		 * SUBROUTINE OPTSTP(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX, +
		 * ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,IPR,MSG)
		 * 
		 * C C UNCONSTRAINED MINIMIZATION STOPPING CRITERIA C
		 * -------------------------------------------- C FIND WHETHER THE
		 * ALGORITHM SHOULD TERMINATE, DUE TO ANY C OF THE FOLLOWING: C 1)
		 * PROBLEM SOLVED WITHIN USER TOLERANCE C 2) CONVERGENCE WITHIN USER
		 * TOLERANCE C 3) ITERATION LIMIT REACHED C 4) DIVERGENCE OR TOO
		 * RESTRICTIVE MAXIMUM STEP (STEPMX) SUSPECTED C C C PARAMETERS C
		 * ---------- C N --> DIMENSION OF PROBLEM C XPLS(N) --> NEW ITERATE
		 * X[K] C FPLS --> FUNCTION VALUE AT NEW ITERATE F(XPLS) C GPLS(N) -->
		 * GRADIENT AT NEW ITERATE, G(XPLS), OR APPROXIMATE C X(N) --> OLD
		 * ITERATE X[K-1] C ITNCNT --> CURRENT ITERATION K C ICSCMX <--> NUMBER
		 * CONSECUTIVE STEPS .GE. STEPMX C [RETAIN VALUE BETWEEN SUCCESSIVE
		 * CALLS] C ITRMCD <-- TERMINATION CODE C GRADTL --> TOLERANCE AT WHICH
		 * RELATIVE GRADIENT CONSIDERED CLOSE C ENOUGH TO ZERO TO TERMINATE
		 * ALGORITHM C STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE
		 * ITERATES C CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM C SX(N) -->
		 * DIAGONAL SCALING MATRIX FOR X C FSCALE --> ESTIMATE OF SCALE OF
		 * OBJECTIVE FUNCTION C ITNLIM --> MAXIMUM NUMBER OF ALLOWABLE
		 * ITERATIONS C IRETCD --> RETURN CODE C MXTAKE --> BOOLEAN FLAG
		 * INDICATING STEP OF MAXIMUM LENGTH USED C IPR --> DEVICE TO WHICH TO
		 * SEND OUTPUT C MSG --> IF MSG INCLUDES A TERM 8, SUPPRESS OUTPUT C C
		 */

		int i;
		double d, rgx, relgrd, rsx, relstp;

		itrmcd[1] = 0;

		// LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X

		if (iretcd[1] == 1) {

			itrmcd[1] = 3;

			if ((msg[1] / 8) % 2 == 0) {

				System.out.print("\n\nOPTSTP    The last global step failed");
				System.out.print(" to locate a point lower than x.\n");
				System.out.print("OPTSTP    Either x is an approximate local");
				System.out.print(" minimum of the function,\n");
				System.out.print("OPTSTP    the function is too nonlinear for");
				System.out.print(" this algorithm, or\n");
				System.out.print("OPTSTP    steptl is too large.\n");

			}

			return;

		} else {

			// FIND DIRECTION IN WHICH RELATIVE GRADIENT MAXIMUM.
			// CHECK WHETHER WITHIN TOLERANCE

			d = Math.max(Math.abs(fpls[1]), fscale[1]);
			rgx = 0.0;

			for (i = 1; i <= n; i++) {

				relgrd = Math.abs(gpls[i])
						* Math.max(Math.abs(xpls[i]), 1.0 / sx[i]) / d;
				rgx = Math.max(rgx, relgrd);

			}

			if (rgx <= gradtl[1]) {

				itrmcd[1] = 1;

				if ((msg[1] / 8) % 2 == 0) {

					System.out
							.print("\n\nOPTSTP    The relative gradient is close");
					System.out.print(" to zero.\n");
					System.out
							.print("OPTSTP    The current iterate is probably");
					System.out.print(" a solution.\n");

				}

				return;

			}

			if (itncnt[1] == 0)
				return;

			// FIND DIRECTION IN WHICH RELATIVE STEPSIZE MAXIMUM
			// CHECK WHETHER WITHIN TOLERANCE.

			rsx = 0.0;

			for (i = 1; i <= n; i++) {

				relstp = Math.abs(xpls[i] - x[i])
						/ Math.max(Math.abs(xpls[i]), 1.0 / sx[i]);

				rsx = Math.max(rsx, relstp);

			}

			if (rsx <= steptl[1]) {

				itrmcd[1] = 2;

				if ((msg[1] / 8) % 2 == 0) {

					System.out
							.print("\n\nOPTSTP    Successive iterates are within");
					System.out.print(" steptl.\n");
					System.out
							.print("OPTSTP    The current iterate is probably");
					System.out.print(" a solution.\n");

				}

				return;

			}

			// CHECK ITERATION LIMIT

			if (itncnt[1] >= itnlim[1]) {

				itrmcd[1] = 4;

				if ((msg[1] / 8) % 2 == 0) {

					System.out
							.print("\n\nOPTSTP    The iteration limit was reached.\n");
					System.out.print("OPTSTP    The algorithm failed.\n");

				}

				return;

			}

			// CHECK NUMBER OF CONSECUTIVE STEPS \ STEPMX

			if (!mxtake[1]) {

				icscmx[1] = 0;

				return;

			}

			if ((msg[1] / 8) % 2 == 0) {

				System.out
						.print("\n\nOPTSTP    Step of maximum length (stepmx)");
				System.out.print(" taken.\n");

			}

			icscmx[1]++;

			if (icscmx[1] < 5)
				return;

			itrmcd[1] = 5;

			if ((msg[1] / 8) % 2 == 0) {

				System.out.print("\n\nOPTSTP    Maximum step size exceeded");
				System.out.print(" five consecutive times.\n");
				System.out.print("OPTSTP    Either the function is unbounded");
				System.out.print(" below,\n");
				System.out
						.print("OPTSTP    becomes asymptotic to a finite value");
				System.out.print(" from above in some direction, or\n");
				System.out.print("OPTSTP    stepmx is too small.\n");

			}

			return;

		}

	}

	/**
	 * 
	 *<p>
	 * The qraux1_f77 method interchanges rows i,i+1 of the upper Hessenberg
	 * matrix r, columns i to n.
	 * 
	 *Translated by Steve Verrill, April 29, 1998.
	 * 
	 *@param n
	 *            The dimension of the matrix
	 *@param r
	 *            Upper Hessenberg matrix
	 *@param i
	 *            Index of row to interchange (i < n)
	 * 
	 * 
	 */

	public static void qraux1_f77(int n, double r[][], int i) {

		/*
		 * 
		 * Here is a copy of the qraux1 FORTRAN documentation:
		 * 
		 * SUBROUTINE QRAUX1(NR,N,R,I) C C PURPOSE C ------- C INTERCHANGE ROWS
		 * I,I+1 OF THE UPPER HESSENBERG MATRIX R, C COLUMNS I TO N C C
		 * PARAMETERS C ---------- C NR --> ROW DIMENSION OF MATRIX C N -->
		 * DIMENSION OF MATRIX C R(N,N) <--> UPPER HESSENBERG MATRIX C I -->
		 * INDEX OF ROW TO INTERCHANGE (I.LT.N) C
		 */

		double tmp;
		int j, ip1;

		ip1 = i + 1;

		for (j = i; j <= n; j++) {

			tmp = r[i][j];
			r[i][j] = r[ip1][j];
			r[ip1][j] = tmp;

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The qraux2_f77 method pre-multiplies r by the Jacobi rotation
	 * j(i,i+1,a,b).
	 * 
	 *Translated by Steve Verrill, April 29, 1998.
	 * 
	 *@param n
	 *            The dimension of the matrix
	 *@param r
	 *            Upper Hessenberg matrix
	 *@param i
	 *            Index of row
	 *@param a
	 *            scalar
	 *@param b
	 *            scalar
	 * 
	 * 
	 */

	public static void qraux2_f77(int n, double r[][], int i, double a, double b) {

		/*
		 * 
		 * Here is a copy of the qraux2 FORTRAN documentation:
		 * 
		 * SUBROUTINE QRAUX2(NR,N,R,I,A,B) C C PURPOSE C ------- C PRE-MULTIPLY
		 * R BY THE JACOBI ROTATION J(I,I+1,A,B) C C PARAMETERS C ---------- C
		 * NR --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF MATRIX C R(N,N)
		 * <--> UPPER HESSENBERG MATRIX C I --> INDEX OF ROW C A --> SCALAR C B
		 * --> SCALAR C
		 */

		double den, c, s, y, z;
		int j, ip1;

		ip1 = i + 1;

		den = Math.sqrt(a * a + b * b);
		c = a / den;
		s = b / den;

		for (j = i; j <= n; j++) {

			y = r[i][j];
			z = r[ip1][j];
			r[i][j] = c * y - s * z;
			r[ip1][j] = s * y + c * z;

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The qrupdt_f77 method finds an orthogonal n by n matrix, Q*, and an upper
	 * triangular n by n matrix, R*, such that (Q*)(R*) = R+U(V+).
	 * 
	 *Translated by Steve Verrill, May 11, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param a
	 *            On input: contains R On output: contains R*
	 *@param u
	 *            Vector
	 *@param v
	 *            Vector
	 * 
	 * 
	 */

	public static void qrupdt_f77(int n, double a[][], double u[], double v[]) {

		/*
		 * 
		 * Here is a copy of the qrupdt FORTRAN documentation:
		 * 
		 * 
		 * SUBROUTINE QRUPDT(NR,N,A,U,V) C C PURPOSE C ------- C FIND AN
		 * ORTHOGONAL (N*N) MATRIX (Q*) AND AN UPPER TRIANGULAR (N*N) C MATRIX
		 * (R*) SUCH THAT (Q*)(R*)=R+U(V+) C C PARAMETERS C ---------- C NR -->
		 * ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C A(N,N) <--> ON
		 * INPUT: CONTAINS R C ON OUTPUT: CONTAINS (R*) C U(N) --> VECTOR C V(N)
		 * --> VECTOR C
		 */

		int k, km1, ii, i, j;
		double t1, t2;

		// DETERMINE LAST NON-ZERO IN U(.)

		k = n;

		while (u[k] == 0 && k > 1) {

			k--;

		}

		// (K-1) JACOBI ROTATIONS TRANSFORM
		// R + U(V+) --> (R*) + (U(1)*E1)(V+)
		// WHICH IS UPPER HESSENBERG

		km1 = k - 1;

		for (ii = 1; ii <= km1; ii++) {

			i = km1 - ii + 1;

			if (u[i] == 0.0) {

				Uncmin_f77.qraux1_f77(n, a, i);
				u[i] = u[i + 1];

			} else {

				Uncmin_f77.qraux2_f77(n, a, i, u[i], -u[i + 1]);
				u[i] = Math.sqrt(u[i] * u[i] + u[i + 1] * u[i + 1]);

			}

		}

		// R <-- R + (U(1)*E1)(V+)

		for (j = 1; j <= n; j++) {

			a[1][j] += u[1] * v[j];

		}

		// (K-1) JACOBI ROTATIONS TRANSFORM UPPER HESSENBERG R
		// TO UPPER TRIANGULAR (R*)

		km1 = k - 1;

		for (i = 1; i <= km1; i++) {

			if (a[i][i] == 0.0) {

				Uncmin_f77.qraux1_f77(n, a, i);

			} else {

				t1 = a[i][i];
				t2 = -a[i + 1][i];

				Uncmin_f77.qraux2_f77(n, a, i, t1, t2);

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The result_f77 method prints information.
	 * 
	 *Translated by Steve Verrill, May 11, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            Estimate of the location of a minimum at iteration k
	 *@param f
	 *            function value at x
	 *@param g
	 *            gradient at x
	 *@param a
	 *            Hessian at x
	 *@param p
	 *            Step taken
	 *@param itncnt
	 *            Iteration number (k)
	 *@param iflg
	 *            Flag controlling the information to print
	 * 
	 */

	public static void result_f77(int n, double x[], double f[], double g[],
			double a[][], double p[], int itncnt[], int iflg) {

		/*
		 * 
		 * Here is a copy of the result FORTRAN documentation:
		 * 
		 * 
		 * SUBROUTINE RESULT(NR,N,X,F,G,A,P,ITNCNT,IFLG,IPR) C C PURPOSE C
		 * ------- C PRINT INFORMATION C C PARAMETERS C ---------- C NR --> ROW
		 * DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> ITERATE
		 * X[K] C F --> FUNCTION VALUE AT X[K] C G(N) --> GRADIENT AT X[K] C
		 * A(N,N) --> HESSIAN AT X[K] C P(N) --> STEP TAKEN C ITNCNT -->
		 * ITERATION NUMBER K C IFLG --> FLAG CONTROLLING INFO TO PRINT C IPR
		 * --> DEVICE TO WHICH TO SEND OUTPUT C
		 */

		int i, j, iii, num5, remain, iii5, iiir;
		int ilow, ihigh;

		num5 = n / 5;
		remain = n % 5;

		// PRINT ITERATION NUMBER

		System.out.print("\n\nRESULT      Iterate k = " + itncnt[1] + "\n");

		if (iflg != 0) {

			// PRINT STEP

			System.out.print("\n\nRESULT      Step\n\n");

			ilow = -4;
			ihigh = 0;

			for (i = 1; i <= num5; i++) {

				ilow += 5;
				ihigh += 5;

				System.out.print(ilow + "--" + ihigh + "     ");

				for (j = 1; j <= 5; j++) {

					System.out.print(p[ilow + j - 1] + "  ");

				}

				System.out.print("\n");

			}

			ilow += 5;
			ihigh = ilow + remain - 1;

			System.out.print(ilow + "--" + ihigh + "     ");

			for (j = 1; j <= remain; j++) {

				System.out.print(p[ilow + j - 1] + "  ");

			}

			System.out.print("\n");

		}

		// PRINT CURRENT ITERATE

		System.out.print("\n\nRESULT      Current x\n\n");

		ilow = -4;
		ihigh = 0;

		for (i = 1; i <= num5; i++) {

			ilow += 5;
			ihigh += 5;

			System.out.print(ilow + "--" + ihigh + "     ");

			for (j = 1; j <= 5; j++) {

				System.out.print(x[ilow + j - 1] + "  ");

			}

			System.out.print("\n");

		}

		ilow += 5;
		ihigh = ilow + remain - 1;

		System.out.print(ilow + "--" + ihigh + "     ");

		for (j = 1; j <= remain; j++) {

			System.out.print(x[ilow + j - 1] + "  ");

		}

		System.out.print("\n");

		// PRINT FUNCTION VALUE

		System.out.print("\n\nRESULT      f_to_minimize at x = " + f[1] + "\n");

		// PRINT GRADIENT

		System.out.print("\n\nRESULT      Gradient at x\n\n");

		ilow = -4;
		ihigh = 0;

		for (i = 1; i <= num5; i++) {

			ilow += 5;
			ihigh += 5;

			System.out.print(ilow + "--" + ihigh + "     ");

			for (j = 1; j <= 5; j++) {

				System.out.print(g[ilow + j - 1] + "  ");

			}

			System.out.print("\n");

		}

		ilow += 5;
		ihigh = ilow + remain - 1;

		System.out.print(ilow + "--" + ihigh + "     ");

		for (j = 1; j <= remain; j++) {

			System.out.print(g[ilow + j - 1] + "  ");

		}

		System.out.print("\n");

		// PRINT HESSIAN FROM ITERATION K

		if (iflg != 0) {

			System.out.print("\n\nRESULT      Hessian at x\n\n");

			for (iii = 1; iii <= n; iii++) {

				iii5 = iii / 5;
				iiir = iii % 5;

				ilow = -4;
				ihigh = 0;

				for (i = 1; i <= iii5; i++) {

					ilow += 5;
					ihigh += 5;

					System.out.print("i = " + iii + ", j = ");
					System.out.print(ilow + "--" + ihigh + "     ");

					for (j = 1; j <= 5; j++) {

						System.out.print(a[iii][ilow + j - 1] + "  ");

					}

					System.out.print("\n");

				}

				ilow += 5;
				ihigh = ilow + iiir - 1;

				System.out.print("i = " + iii + ", j = ");
				System.out.print(ilow + "--" + ihigh + "     ");

				for (j = 1; j <= iiir; j++) {

					System.out.print(a[iii][ilow + j - 1] + "  ");

				}

				System.out.print("\n");

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The sclmul_f77 method multiplies a vector by a scalar.
	 * 
	 *Translated by Steve Verrill, May 8, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param s
	 *            The scalar
	 *@param v
	 *            Operand vector
	 *@param z
	 *            Result vector
	 * 
	 * 
	 */

	public static void sclmul_f77(int n, double s, double v[], double z[]) {

		/*
		 * 
		 * Here is a copy of the sclmul FORTRAN documentation:
		 * 
		 * SUBROUTINE SCLMUL(N,S,V,Z) C C PURPOSE C ------- C MULTIPLY VECTOR BY
		 * SCALAR C RESULT VECTOR MAY BE OPERAND VECTOR C C PARAMETERS C
		 * ---------- C N --> DIMENSION OF VECTORS C S --> SCALAR C V(N) -->
		 * OPERAND VECTOR C Z(N) <-- RESULT VECTOR
		 */

		int i;

		for (i = 1; i <= n; i++) {

			z[i] = s * v[i];

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The secfac_f77 method updates the Hessian by the BFGS factored technique.
	 * 
	 *Translated by Steve Verrill, May 14, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            Old iterate
	 *@param g
	 *            Gradient or approximation at the old iterate
	 *@param a
	 *            On entry: Cholesky decomposition of Hessian in lower triangle
	 *            and diagonal On exit: Updated Cholesky decomposition of
	 *            Hessian in lower triangle and diagonal
	 *@param xpls
	 *            New iterate
	 *@param gpls
	 *            Gradient or approximation at the new iterate
	 *@param epsm
	 *            Machine epsilon
	 *@param itncnt
	 *            Iteration count
	 *@param rnf
	 *            Relative noise in optimization function f_to_minimize
	 *@param iagflg
	 *            1 if an analytic gradient is supplied, 0 otherwise
	 *@param noupdt
	 *            Boolean: no update yet (retain value between successive calls)
	 *@param s
	 *            Workspace
	 *@param y
	 *            Workspace
	 *@param u
	 *            Workspace
	 *@param w
	 *            Workspace
	 * 
	 * 
	 */

	public static void secfac_f77(int n, double x[], double g[], double a[][],
			double xpls[], double gpls[], double epsm, int itncnt[],
			double rnf, int iagflg[], boolean noupdt[], double s[], double y[],
			double u[], double w[]) {

		/*
		 * 
		 * Here is a copy of the secfac FORTRAN documentation:
		 * 
		 * SUBROUTINE SECFAC(NR,N,X,G,A,XPLS,GPLS,EPSM,ITNCNT,RNF, +
		 * IAGFLG,NOUPDT,S,Y,U,W) C C PURPOSE C ------- C UPDATE HESSIAN BY THE
		 * BFGS FACTORED METHOD C C PARAMETERS C ---------- C NR --> ROW
		 * DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> OLD
		 * ITERATE, X[K-1] C G(N) --> GRADIENT OR APPROXIMATE AT OLD ITERATE C
		 * A(N,N) <--> ON ENTRY: CHOLESKY DECOMPOSITION OF HESSIAN IN C LOWER
		 * PART AND DIAGONAL. C ON EXIT: UPDATED CHOLESKY DECOMPOSITION OF
		 * HESSIAN C IN LOWER TRIANGULAR PART AND DIAGONAL C XPLS(N) --> NEW
		 * ITERATE, X[K] C GPLS(N) --> GRADIENT OR APPROXIMATE AT NEW ITERATE C
		 * EPSM --> MACHINE EPSILON C ITNCNT --> ITERATION COUNT C RNF -->
		 * RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN C IAGFLG --> =1 IF
		 * ANALYTIC GRADIENT SUPPLIED, =0 ITHERWISE C NOUPDT <--> BOOLEAN: NO
		 * UPDATE YET C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C S(N) -->
		 * WORKSPACE C Y(N) --> WORKSPACE C U(N) --> WORKSPACE C W(N) -->
		 * WORKSPACE C
		 */

		boolean skpupd;
		int i, j, im1;
		double den1, snorm2, ynrm2, den2, alp, reltol;

		if (itncnt[1] == 1)
			noupdt[1] = true;

		for (i = 1; i <= n; i++) {

			s[i] = xpls[i] - x[i];
			y[i] = gpls[i] - g[i];

		}

		den1 = Blas_f77.ddot_f77(n, s, 1, y, 1);
		snorm2 = Blas_f77.dnrm2_f77(n, s, 1);
		ynrm2 = Blas_f77.dnrm2_f77(n, y, 1);

		if (den1 >= Math.sqrt(epsm) * snorm2 * ynrm2) {

			Uncmin_f77.mvmltu_f77(n, a, s, u);
			den2 = Blas_f77.ddot_f77(n, u, 1, u, 1);

			// L <-- SQRT(DEN1/DEN2)*L

			alp = Math.sqrt(den1 / den2);

			if (noupdt[1]) {

				for (j = 1; j <= n; j++) {

					u[j] *= alp;

					for (i = j; i <= n; i++) {

						a[i][j] *= alp;

					}

				}

				noupdt[1] = false;
				den2 = den1;
				alp = 1.0;

			}

			skpupd = true;

			// W = L(L+)S = HS

			Uncmin_f77.mvmltl_f77(n, a, u, w);

			i = 1;

			if (iagflg[1] == 0) {

				reltol = Math.sqrt(rnf);

			} else {

				reltol = rnf;

			}

			while (i <= n && skpupd) {

				if (Math.abs(y[i] - w[i]) >= reltol
						* Math.max(Math.abs(g[i]), Math.abs(gpls[i]))) {

					skpupd = false;

				} else {

					i++;

				}

			}

			if (!skpupd) {

				// W=Y-ALP*L(L+)S

				for (i = 1; i <= n; i++) {

					w[i] = y[i] - alp * w[i];

				}

				// ALP=1/SQRT(DEN1*DEN2)

				alp /= den1;

				// U=(L+)/SQRT(DEN1*DEN2) = (L+)S/SQRT((Y+)S * (S+)L(L+)S)

				for (i = 1; i <= n; i++) {

					u[i] *= alp;

				}

				// COPY L INTO UPPER TRIANGULAR PART. ZERO L.

				for (i = 2; i <= n; i++) {

					im1 = i - 1;

					for (j = 1; j <= im1; j++) {

						a[j][i] = a[i][j];
						a[i][j] = 0.0;

					}

				}

				// FIND Q, (L+) SUCH THAT Q(L+) = (L+) + U(W+)

				Uncmin_f77.qrupdt_f77(n, a, u, w);

				// UPPER TRIANGULAR PART AND DIAGONAL OF A NOW CONTAIN UPDATED
				// CHOLESKY DECOMPOSITION OF HESSIAN. COPY BACK TO LOWER
				// TRIANGULAR PART.

				for (i = 2; i <= n; i++) {

					im1 = i - 1;

					for (j = 1; j <= im1; j++) {

						a[i][j] = a[j][i];

					}

				}

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The secunf_f77 method updates the Hessian by the BFGS unfactored
	 * approach.
	 * 
	 *Translated by Steve Verrill, May 8, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            The old iterate
	 *@param g
	 *            The gradient or an approximation at the old iterate
	 *@param a
	 *            On entry: Approximate Hessian at the old iterate in the upper
	 *            triangular part (and udiag) On exit: Updated approximate
	 *            Hessian at the new iterate in the lower triangular part and
	 *            diagonal
	 *@param udiag
	 *            On entry: Diagonal of Hessian
	 *@param xpls
	 *            New iterate
	 *@param gpls
	 *            Gradient or approximation at the new iterate
	 *@param epsm
	 *            Machine epsilon
	 *@param itncnt
	 *            Iteration count
	 *@param rnf
	 *            Relative noise in the optimization function, f_to_minimize
	 *@param iagflg
	 *            = 1 if an analytic gradient is supplied, = 0 otherwise
	 *@param noupdt
	 *            Boolean: no update yet (retain value between calls)
	 *@param s
	 *            workspace
	 *@param y
	 *            workspace
	 *@param t
	 *            workspace
	 * 
	 * 
	 */

	public static void secunf_f77(int n, double x[], double g[], double a[][],
			double udiag[], double xpls[], double gpls[], double epsm,
			int itncnt[], double rnf, int iagflg[], boolean noupdt[],
			double s[], double y[], double t[]) {

		/*
		 * 
		 * Here is a copy of the secunf FORTRAN documentation:
		 * 
		 * SUBROUTINE SECUNF(NR,N,X,G,A,UDIAG,XPLS,GPLS,EPSM,ITNCNT, +
		 * RNF,IAGFLG,NOUPDT,S,Y,T) C C PURPOSE C ------- C UPDATE HESSIAN BY
		 * THE BFGS UNFACTORED METHOD C C PARAMETERS C ---------- C NR --> ROW
		 * DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) --> OLD
		 * ITERATE, X[K-1] C G(N) --> GRADIENT OR APPROXIMATE AT OLD ITERATE C
		 * A(N,N) <--> ON ENTRY: APPROXIMATE HESSIAN AT OLD ITERATE C IN UPPER
		 * TRIANGULAR PART (AND UDIAG) C ON EXIT: UPDATED APPROX HESSIAN AT NEW
		 * ITERATE C IN LOWER TRIANGULAR PART AND DIAGONAL C [LOWER TRIANGULAR
		 * PART OF SYMMETRIC MATRIX] C UDIAG --> ON ENTRY: DIAGONAL OF HESSIAN C
		 * XPLS(N) --> NEW ITERATE, X[K] C GPLS(N) --> GRADIENT OR APPROXIMATE
		 * AT NEW ITERATE C EPSM --> MACHINE EPSILON C ITNCNT --> ITERATION
		 * COUNT C RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN C IAGFLG
		 * --> =1 IF ANALYTIC GRADIENT SUPPLIED, =0 OTHERWISE C NOUPDT <-->
		 * BOOLEAN: NO UPDATE YET C [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C
		 * S(N) --> WORKSPACE C Y(N) --> WORKSPACE C T(N) --> WORKSPACE C
		 */

		double den1, snorm2, ynrm2, den2, gam, tol;
		int i, j;
		boolean skpupd;

		// COPY HESSIAN IN UPPER TRIANGULAR PART AND UDIAG TO
		// LOWER TRIANGULAR PART AND DIAGONAL

		for (j = 1; j <= n; j++) {

			a[j][j] = udiag[j];

			for (i = j + 1; i <= n; i++) {

				a[i][j] = a[j][i];

			}

		}

		if (itncnt[1] == 1)
			noupdt[1] = true;

		for (i = 1; i <= n; i++) {

			s[i] = xpls[i] - x[i];
			y[i] = gpls[i] - g[i];

		}

		den1 = Blas_f77.ddot_f77(n, s, 1, y, 1);

		snorm2 = Blas_f77.dnrm2_f77(n, s, 1);

		ynrm2 = Blas_f77.dnrm2_f77(n, y, 1);

		if (den1 >= Math.sqrt(epsm) * snorm2 * ynrm2) {

			Uncmin_f77.mvmlts_f77(n, a, s, t);

			den2 = Blas_f77.ddot_f77(n, s, 1, t, 1);

			if (noupdt[1]) {

				// H <-- [(S+)Y/(S+)HS]H

				gam = den1 / den2;

				den2 = gam * den2;

				for (j = 1; j <= n; j++) {

					t[j] *= gam;

					for (i = j; i <= n; i++) {

						a[i][j] *= gam;

					}

				}

				noupdt[1] = false;

			}

			skpupd = true;

			// CHECK UPDATE CONDITION ON ROW I

			for (i = 1; i <= n; i++) {

				tol = rnf * Math.max(Math.abs(g[i]), Math.abs(gpls[i]));
				if (iagflg[1] == 0)
					tol /= Math.sqrt(rnf);

				if (Math.abs(y[i] - t[i]) >= tol) {

					skpupd = false;
					break;

				}

			}

			if (!skpupd) {

				// BFGS UPDATE

				for (j = 1; j <= n; j++) {

					for (i = j; i <= n; i++) {

						a[i][j] += y[i] * y[j] / den1 - t[i] * t[j] / den2;

					}

				}

			}

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The sndofd_f77 method finds second order forward finite difference
	 * approximations to the Hessian. For optimization use this method to
	 * estimate the Hessian of the optimization function if no analytical user
	 * function has been supplied for either the gradient or the Hessian, and
	 * the optimization function is inexpensive to evaluate.
	 * 
	 *Translated by Steve Verrill, May 8, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param xpls
	 *            New iterate
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param fpls
	 *            Function value at the new iterate
	 *@param a
	 *            "FINITE DIFFERENCE APPROXIMATION TO HESSIAN. ONLY LOWER
	 *            TRIANGULAR MATRIX AND DIAGONAL ARE RETURNED"
	 *@param sx
	 *            Scaling vector for x
	 *@param rnoise
	 *            Relative noise in the function to be minimized
	 *@param stepsz
	 *            Workspace (stepsize in i-th component direction)
	 *@param anbr
	 *            Workspace (neighbor in i-th direction)
	 * 
	 */

	public static void sndofd_f77(int n, double xpls[],
			Uncmin_methods minclass, double fpls[], double a[][], double sx[],
			double rnoise, double stepsz[], double anbr[]) {

		/*
		 * 
		 * Here is a copy of the sndofd FORTRAN documentation.
		 * 
		 * SUBROUTINE SNDOFD(NR,N,XPLS,FCN,FPLS,A,SX,RNOISE,STEPSZ,ANBR)
		 * 
		 * C PURPOSE C ------- C FIND SECOND ORDER FORWARD FINITE DIFFERENCE
		 * APPROXIMATION "A" C TO THE SECOND DERIVATIVE (HESSIAN) OF THE
		 * FUNCTION DEFINED BY THE SUBP C "FCN" EVALUATED AT THE NEW ITERATE
		 * "XPLS" C C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE C 1) THE
		 * SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION C IF NO
		 * ANALYTICAL USER FUNCTION HAS BEEN SUPPLIED FOR EITHER C THE GRADIENT
		 * OR THE HESSIAN AND IF THE OPTIMIZATION FUNCTION C "FCN" IS
		 * INEXPENSIVE TO EVALUATE. C C PARAMETERS C ---------- C NR --> ROW
		 * DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C XPLS(N) --> NEW
		 * ITERATE: X[K] C FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION C
		 * FPLS --> FUNCTION VALUE AT NEW ITERATE, F(XPLS) C A(N,N) <-- FINITE
		 * DIFFERENCE APPROXIMATION TO HESSIAN C ONLY LOWER TRIANGULAR MATRIX
		 * AND DIAGONAL C ARE RETURNED C SX(N) --> DIAGONAL SCALING MATRIX FOR X
		 * C RNOISE --> RELATIVE NOISE IN FNAME [F(X)] C STEPSZ(N) --> WORKSPACE
		 * (STEPSIZE IN I-TH COMPONENT DIRECTION) C ANBR(N) --> WORKSPACE
		 * (NEIGHBOR IN I-TH DIRECTION) C C
		 */

		double xmult, xtmpi, xtmpj, fhat;
		int i, j;

		// FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION
		// OF I-TH UNIT VECTOR

		xmult = Math.pow(rnoise, 1.0 / 3.0);

		for (i = 1; i <= n; i++) {

			stepsz[i] = xmult * Math.max(Math.abs(xpls[i]), 1.0 / sx[i]);
			xtmpi = xpls[i];
			xpls[i] = xtmpi + stepsz[i];

			anbr[i] = minclass.f_to_minimize(xpls);
			xpls[i] = xtmpi;

		}

		// CALCULATE COLUMN I OF A

		for (i = 1; i <= n; i++) {

			xtmpi = xpls[i];
			xpls[i] = xtmpi + 2.0 * stepsz[i];
			fhat = minclass.f_to_minimize(xpls);
			a[i][i] = ((fpls[1] - anbr[i]) + (fhat - anbr[i]))
					/ (stepsz[i] * stepsz[i]);

			// CALCULATE SUB-DIAGONAL ELEMENTS OF COLUMN

			if (i != n) {

				xpls[i] = xtmpi + stepsz[i];

				for (j = i + 1; j <= n; j++) {

					xtmpj = xpls[j];
					xpls[j] = xtmpj + stepsz[j];
					fhat = minclass.f_to_minimize(xpls);
					a[j][i] = ((fpls[1] - anbr[i]) + (fhat - anbr[j]))
							/ (stepsz[i] * stepsz[j]);
					xpls[j] = xtmpj;

				}

			}

			xpls[i] = xtmpi;

		}

		return;

	}

	/**
	 * 
	 *<p>
	 * The tregup_f77 method decides whether to accept xpls = x + sc as the next
	 * iterate and update the trust region dlt.
	 * 
	 *Translated by Steve Verrill, May 11, 1998.
	 * 
	 *@param n
	 *            The dimension of the problem
	 *@param x
	 *            Old iterate
	 *@param f
	 *            Function value at old iterate
	 *@param g
	 *            Gradient or approximation at old iterate
	 *@param a
	 *            Cholesky decomposition of Hessian in lower triangular part and
	 *            diagonal. Hessian or approximation in upper triangular part.
	 *@param minclass
	 *            A class that implements the Uncmin_methods interface (see the
	 *            definition in Uncmin_methods.java). See UncminTest_f77.java
	 *            for an example of such a class. The class must define: 1.) a
	 *            method, f_to_minimize, to minimize. f_to_minimize must have
	 *            the form
	 * 
	 *            public static double f_to_minimize(double x[])
	 * 
	 *            where x is the vector of arguments to the function and the
	 *            return value is the value of the function evaluated at x. 2.)
	 *            a method, gradient, that has the form
	 * 
	 *            public static void gradient(double x[], double g[])
	 * 
	 *            where g is the gradient of f evaluated at x. This method will
	 *            have an empty body if the user does not wish to provide an
	 *            analytic estimate of the gradient. 3.) a method, hessian, that
	 *            has the form
	 * 
	 *            public static void hessian(double x[], double h[][]) where h
	 *            is the Hessian of f evaluated at x. This method will have an
	 *            empty body if the user does not wish to provide an analytic
	 *            estimate of the Hessian. If the user wants Uncmin to check the
	 *            Hessian, then the hessian method should only fill the lower
	 *            triangle (and diagonal) of h.
	 *@param sc
	 *            Current step
	 *@param sx
	 *            Scaling vector for x
	 *@param nwtake
	 *            Boolean, = true if Newton step taken
	 *@param stepmx
	 *            Maximum allowable step size
	 *@param steptl
	 *            Relative step size at which successive iterates are considered
	 *            close enough to terminate the algorithm
	 *@param dlt
	 *            Trust region radius
	 *@param iretcd
	 *            Return code = 0 xpls accepted as next iterate, dlt is the
	 *            trust region radius for the next iteration = 1 xpls
	 *            unsatisfactory but accepted as next iterate because xpls - x
	 *            is less than the smallest allowable step length = 2 f(xpls)
	 *            too large. Continue current iteration with new reduced dlt. =
	 *            3 f(xpls) sufficiently small, but quadratic model predicts
	 *            f(xpls) sufficiently well to continue current iteration with
	 *            new doubled dlt.
	 *@param xplsp
	 *            Workspace (value needs to be retained between successive calls
	 *            of k-th global step)
	 *@param fplsp
	 *            Retain between successive calls
	 *@param xpls
	 *            New iterate
	 *@param fpls
	 *            Function value at new iterate
	 *@param mxtake
	 *            Boolean flag indicating step of maximum length used
	 *@param method
	 *            Algorithm to use to solve minimization problem = 1 Line search
	 *            = 2 Double dogleg = 3 More-Hebdon
	 *@param udiag
	 *            Diagonal of Hessian in a
	 * 
	 */

	public static void tregup_f77(int n, double x[], double f[], double g[],
			double a[][], Uncmin_methods minclass, double sc[], double sx[],
			boolean nwtake[], double stepmx[], double steptl[], double dlt[],
			int iretcd[], double xplsp[], double fplsp[], double xpls[],
			double fpls[], boolean mxtake[], int method, double udiag[]) {

		/*
		 * 
		 * Here is a copy of the tregup FORTRAN documentation.
		 * 
		 * 
		 * SUBROUTINE TREGUP(NR,N,X,F,G,A,FCN,SC,SX,NWTAKE,STEPMX,STEPTL, +
		 * DLT,IRETCD,XPLSP,FPLSP,XPLS,FPLS,MXTAKE,IPR,METHOD,UDIAG) C C PURPOSE
		 * C ------- C DECIDE WHETHER TO ACCEPT XPLS=X+SC AS THE NEXT ITERATE
		 * AND UPDATE THE C TRUST REGION DLT. C C PARAMETERS C ---------- C NR
		 * --> ROW DIMENSION OF MATRIX C N --> DIMENSION OF PROBLEM C X(N) -->
		 * OLD ITERATE X[K-1] C F --> FUNCTION VALUE AT OLD ITERATE, F(X) C G(N)
		 * --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE C A(N,N) -->
		 * CHOLESKY DECOMPOSITION OF HESSIAN IN C LOWER TRIANGULAR PART AND
		 * DIAGONAL. C HESSIAN OR APPROX IN UPPER TRIANGULAR PART C FCN --> NAME
		 * OF SUBROUTINE TO EVALUATE FUNCTION C SC(N) --> CURRENT STEP C SX(N)
		 * --> DIAGONAL SCALING MATRIX FOR X C NWTAKE --> BOOLEAN, =.TRUE. IF
		 * NEWTON STEP TAKEN C STEPMX --> MAXIMUM ALLOWABLE STEP SIZE C STEPTL
		 * --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES C CONSIDERED
		 * CLOSE ENOUGH TO TERMINATE ALGORITHM C DLT <--> TRUST REGION RADIUS C
		 * IRETCD <--> RETURN CODE C =0 XPLS ACCEPTED AS NEXT ITERATE; C DLT
		 * TRUST REGION FOR NEXT ITERATION. C =1 XPLS UNSATISFACTORY BUT
		 * ACCEPTED AS NEXT ITERATE C BECAUSE XPLS-X .LT. SMALLEST ALLOWABLE C
		 * STEP LENGTH. C =2 F(XPLS) TOO LARGE. CONTINUE CURRENT ITERATION C
		 * WITH NEW REDUCED DLT. C =3 F(XPLS) SUFFICIENTLY SMALL, BUT QUADRATIC
		 * MODEL C PREDICTS F(XPLS) SUFFICIENTLY WELL TO CONTINUE C CURRENT
		 * ITERATION WITH NEW DOUBLED DLT. C XPLSP(N) <--> WORKSPACE [VALUE
		 * NEEDS TO BE RETAINED BETWEEN C SUCCESIVE CALLS OF K-TH GLOBAL STEP] C
		 * FPLSP <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] C XPLS(N) <-- NEW
		 * ITERATE X[K] C FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS) C
		 * MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED C IPR
		 * --> DEVICE TO WHICH TO SEND OUTPUT C METHOD --> ALGORITHM TO USE TO
		 * SOLVE MINIMIZATION PROBLEM C =1 LINE SEARCH C =2 DOUBLE DOGLEG C =3
		 * MORE-HEBDON C UDIAG(N) --> DIAGONAL OF HESSIAN IN A(.,.) C
		 */

		int i, j;
		double rln, temp, dltf, slp, dltmp, dltfp;

		mxtake[1] = false;

		for (i = 1; i <= n; i++) {

			xpls[i] = x[i] + sc[i];

		}

		fpls[1] = minclass.f_to_minimize(xpls);

		dltf = fpls[1] - f[1];

		slp = Blas_f77.ddot_f77(n, g, 1, sc, 1);

		if (iretcd[1] == 4)
			fplsp[1] = 0.0;

		if ((iretcd[1] == 3) && ((fpls[1] >= fplsp[1]) || (dltf > .0001 * slp))) {

			// RESET XPLS TO XPLSP AND TERMINATE GLOBAL STEP

			iretcd[1] = 0;

			for (i = 1; i <= n; i++) {

				xpls[i] = xplsp[i];

			}

			fpls[1] = fplsp[1];
			dlt[1] *= .5;

		} else {

			// FPLS TOO LARGE

			if (dltf > .0001 * slp) {

				rln = 0.0;

				for (i = 1; i <= n; i++) {

					rln = Math.max(rln, Math.abs(sc[i])
							/ Math.max(Math.abs(xpls[i]), 1.0 / sx[i]));

				}

				if (rln < steptl[1]) {

					// CANNOT FIND SATISFACTORY XPLS SUFFICIENTLY DISTINCT FROM
					// X

					iretcd[1] = 1;

				} else {

					// REDUCE TRUST REGION AND CONTINUE GLOBAL STEP

					iretcd[1] = 2;

					dltmp = -slp * dlt[1] / (2.0 * (dltf - slp));

					if (dltmp < .1 * dlt[1]) {

						dlt[1] *= .1;

					} else {

						dlt[1] = dltmp;

					}

				}

			} else {

				// FPLS SUFFICIENTLY SMALL

				dltfp = 0.0;

				if (method == 2) {

					for (i = 1; i <= n; i++) {

						temp = 0.0;

						for (j = i; j <= n; j++) {

							temp += a[j][i] * sc[j];

						}

						dltfp += temp * temp;

					}

				} else {

					for (i = 1; i <= n; i++) {

						dltfp += udiag[i] * sc[i] * sc[i];
						temp = 0.0;

						for (j = i + 1; j <= n; j++) {

							temp += a[i][j] * sc[i] * sc[j];

						}

						dltfp += 2.0 * temp;

					}

				}

				dltfp = slp + dltfp / 2.0;

				if ((iretcd[1] != 2)
						&& (Math.abs(dltfp - dltf) <= .1 * Math.abs(dltf))
						&& (!nwtake[1]) && (dlt[1] <= .99 * stepmx[1])) {

					// DOUBLE TRUST REGION AND CONTINUE GLOBAL STEP

					iretcd[1] = 3;

					for (i = 1; i <= n; i++) {

						xplsp[i] = xpls[i];

					}

					fplsp[1] = fpls[1];
					dlt[1] = Math.min(2.0 * dlt[1], stepmx[1]);

				} else {

					// ACCEPT XPLS AS NEXT ITERATE. CHOOSE NEW TRUST REGION.

					iretcd[1] = 0;

					if (dlt[1] > .99 * stepmx[1])
						mxtake[1] = true;

					if (dltf >= .1 * dltfp) {

						// DECREASE TRUST REGION FOR NEXT ITERATION

						dlt[1] *= .5;

					} else {

						// CHECK WHETHER TO INCREASE TRUST REGION FOR NEXT
						// ITERATION

						if (dltf <= .75 * dltfp)
							dlt[1] = Math.min(2.0 * dlt[1], stepmx[1]);

					}

				}

			}

		}

		return;

	}

}
